!> \brief Routines for initialization, allocation, changes to geometry

MODULE INIT

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
USE OUTPUT_DATA
USE TRAN
USE MEMORY_FUNCTIONS, ONLY : CHKMEMERR
USE DEVICE_VARIABLES

IMPLICIT NONE (TYPE,EXTERNAL)

PRIVATE

PUBLIC INITIALIZE_MESH_VARIABLES_1,INITIALIZE_MESH_VARIABLES_2,INITIALIZE_MESH_VARIABLES_3,INITIALIZE_GLOBAL_VARIABLES, &
       OPEN_AND_CLOSE,INITIAL_NOISE,UVW_INIT,TMP_INIT,SPEC_INIT,INITIALIZE_DEVICES,INITIALIZE_PROFILES,REASSIGN_WALL_CELLS,&
       ADJUST_HT3D_WALL_CELLS,INITIALIZE_HT3D_WALL_CELLS,FIND_WALL_BACK_INDICES

CONTAINS


!> \brief Allocate the bulk of arrays used throughout the simulation
!> \param DT Time step (s)
!> \param NM Mesh number

SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)

USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_SPECIFIC_GAS_CONSTANT,GET_SPECIFIC_HEAT,LES_FILTER_WIDTH_FUNCTION,&
                              COMPUTE_WIND_COMPONENTS
USE RADCONS, ONLY: UIIDIM
USE CONTROL_VARIABLES
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER :: N,I,J,K,IW,IC,SURF_INDEX,IOR,IERR,IZERO,II,JJ,KK,OBST_INDEX,N_EXTERNAL_CELLS,NS,N_LOWER_SCALARS
REAL(EB), INTENT(IN) :: DT
INTEGER, INTENT(IN) :: NM
REAL(EB) :: MU_N,CS,DELTA,INTEGRAL,TEMP,ZSW
REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: ZZ_GET,VF
INTEGER, POINTER :: IBP1, JBP1, KBP1,IBAR, JBAR, KBAR
REAL(EB),POINTER :: XS,XF,YS,YF,ZS,ZF
TYPE (INITIALIZATION_TYPE), POINTER :: IN
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (WALL_TYPE), POINTER :: WC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (MESH_TYPE), POINTER :: M
TYPE (RAMPS_TYPE), POINTER :: RP
TYPE (MULTIPLIER_TYPE), POINTER :: MR

IERR = 0
M => MESHES(NM)
IBP1 =>M%IBP1
JBP1 =>M%JBP1
KBP1 =>M%KBP1
IBAR =>M%IBAR
JBAR =>M%JBAR
KBAR =>M%KBAR
XS=>M%XS
YS=>M%YS
ZS=>M%ZS
XF=>M%XF
YF=>M%YF
ZF=>M%ZF

ALLOCATE(M%EXTERNAL_WALL(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','EXTERNAL_WALL',IZERO)

ALLOCATE(M%RHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO',IZERO)
ALLOCATE(M%RHOS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHOS',IZERO)
M%RHOS = RHOA
ALLOCATE(M%TMP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP',IZERO)
ALLOCATE(M%U(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','U',IZERO)
ALLOCATE(M%V(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','V',IZERO)
ALLOCATE(M%W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','W',IZERO)
ALLOCATE(M%US(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','US',IZERO)
ALLOCATE(M%VS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','VS',IZERO)
ALLOCATE(M%WS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','WS',IZERO)
ALLOCATE(M%FVX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX',IZERO)
ALLOCATE(M%FVY(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY',IZERO)
ALLOCATE(M%FVZ(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ',IZERO)
ALLOCATE(M%FVX_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVX_B',IZERO) ; M%FVX_B=0._EB
ALLOCATE(M%FVY_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVY_B',IZERO) ; M%FVY_B=0._EB
ALLOCATE(M%FVZ_B(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','FVZ_B',IZERO) ; M%FVZ_B=0._EB
IF (PARTICLE_DRAG) THEN
   ALLOCATE(M%FVX_D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','FVX_D',IZERO) ; M%FVX_D=0._EB
   ALLOCATE(M%FVY_D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','FVY_D',IZERO) ; M%FVY_D=0._EB
   ALLOCATE(M%FVZ_D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','FVZ_D',IZERO) ; M%FVZ_D=0._EB
ENDIF
ALLOCATE(M%H(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','H',IZERO)
ALLOCATE(M%HS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','HS',IZERO)
ALLOCATE(M%KRES(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','KRES',IZERO)
ALLOCATE(M%DDDT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DDDT',IZERO)
ALLOCATE(M%D(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','D',IZERO)
ALLOCATE(M%DS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','DS',IZERO)
ALLOCATE(M%MU(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU',IZERO)
ALLOCATE(M%MU_DNS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MU_DNS',IZERO); M%MU_DNS = 0._EB
IF (CHECK_VN) THEN
   ALLOCATE(M%D_Z_MAX(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_Z_MAX',IZERO)
   M%D_Z_MAX=0._EB
ENDIF
ALLOCATE(M%STRAIN_RATE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','STRAIN_RATE',IZERO)
M%STRAIN_RATE = 0._EB
SELECT CASE(TURB_MODEL)
   CASE (CONSMAG,DYNSMAG)
      ALLOCATE(M%CSD2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
      CALL ChkMemErr('INIT','CS',IZERO)
END SELECT
IF (OUTPUT_CHEM_IT) THEN
   ALLOCATE(M%CHEM_SUBIT(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','CHEM_SUBIT',IZERO)
   M%CHEM_SUBIT = 0._EB
ENDIF

ALLOCATE(M%CHEM_ACTIVE_CELLS(IBP1*JBP1*KBP1,3),STAT=IZERO)
CALL ChkMemErr('INIT','CHEM_ACTIVE_CELLS',IZERO)
M%CHEM_ACTIVE_CELLS=-1

ALLOCATE(M%Q(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','Q',IZERO)

ALLOCATE(M%MIX_TIME(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','MIX_TIME',IZERO)
M%MIX_TIME=DT

! Background pressure, temperature, density as a function of height (Z coordinate)

ALLOCATE(M%PRESSURE_ZONE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','PRESSURE_ZONE',IZERO) ; M%PRESSURE_ZONE = -1

ALLOCATE(M%P_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','P_0',IZERO)
ALLOCATE(M%TMP_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','TMP_0',IZERO)
ALLOCATE(M%RHO_0(0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RHO_0',IZERO)

! Allocate species arrays

ALLOCATE( M%ZZ(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','ZZ',IZERO)
M%ZZ = 0._EB
ALLOCATE(M%ZZS(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','ZZS',IZERO)
M%ZZS = 0._EB
ALLOCATE(M%DEL_RHO_D_DEL_Z(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','DEL_RHO_D_DEL_Z',IZERO)
M%DEL_RHO_D_DEL_Z = 0._EB
IF (REAC_SOURCE_CHECK) THEN
   ALLOCATE( M%REAC_SOURCE_TERM(0:IBP1,0:JBP1,0:KBP1,N_TRACKED_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','REAC_SOURCE_TERM',IZERO)
   M%REAC_SOURCE_TERM = 0._EB
   ALLOCATE( M%Q_REAC(0:IBP1,0:JBP1,0:KBP1,N_REACTIONS),STAT=IZERO)
   CALL ChkMemErr('INIT','Q_REAC',IZERO)
   M%Q_REAC = 0._EB
ENDIF

ALLOCATE(M%RSUM(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
CALL ChkMemErr('INIT','RSUM',IZERO)
M%RSUM = RSUM0

! Allocate scalar face values

! Required for cell face density correction for multicomponent mixtures
IF (TEST_FLUX_LIMITER_FACE_CORRECTION .AND. N_TRACKED_SPECIES>2) THEN
   N_LOWER_SCALARS=0
ELSE
   N_LOWER_SCALARS=1
ENDIF

ALLOCATE( M%FX(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FX',IZERO)
ALLOCATE( M%FY(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FY',IZERO)
ALLOCATE( M%FZ(0:IBP1,0:JBP1,0:KBP1,N_LOWER_SCALARS:N_TOTAL_SCALARS),STAT=IZERO)
CALL ChkMemErr('INIT','FZ',IZERO)
M%FX = 0._EB
M%FY = 0._EB
M%FZ = 0._EB

! Allocate storage for scalar total fluxes

IF (STORE_SPECIES_FLUX) THEN
   ALLOCATE( M%ADV_FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FX',IZERO)
   ALLOCATE( M%ADV_FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FY',IZERO)
   ALLOCATE( M%ADV_FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','ADV_FZ',IZERO)
   ALLOCATE( M%DIF_FX(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FX',IZERO)
   ALLOCATE( M%DIF_FY(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FY',IZERO)
   ALLOCATE( M%DIF_FZ(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FZ',IZERO)
   ALLOCATE( M%DIF_FXS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FX',IZERO)
   ALLOCATE( M%DIF_FYS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FY',IZERO)
   ALLOCATE( M%DIF_FZS(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO)
   CALL ChkMemErr('INIT','DIF_FZ',IZERO)
   M%ADV_FX = 0._EB
   M%ADV_FY = 0._EB
   M%ADV_FZ = 0._EB
   M%DIF_FX = 0._EB
   M%DIF_FY = 0._EB
   M%DIF_FZ = 0._EB
   M%DIF_FXS = 0._EB
   M%DIF_FYS = 0._EB
   M%DIF_FZS = 0._EB
ENDIF

! Allocate array to store pressure Poisson residual for output

IF (STORE_PRESSURE_POISSON_RESIDUAL) THEN
   ALLOCATE(M%PP_RESIDUAL(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','PP_RESIDUAL',IZERO)
   M%PP_RESIDUAL = 0._EB
ENDIF

! Allocate array to store cut-cell divergence if needed

IF (STORE_CUTCELL_DIVERGENCE) THEN
   ALLOCATE(M%CCVELDIV(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','CCVELDIV',IZERO)
   M%CCVELDIV = 0._EB
ENDIF
IF (STORE_CARTESIAN_DIVERGENCE) THEN
   ALLOCATE(M%CARTVELDIV(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','CARTVELDIV',IZERO)
   M%CARTVELDIV = 0._EB
ENDIF

! Allocate water mass arrays if sprinklers are present

IF (N_LP_ARRAY_INDICES>0) THEN
   ALLOCATE(M%QR_W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','QR_W',IZERO)
   M%QR_W = 0._EB
   ALLOCATE(M%AVG_DROP_DEN(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_DEN',IZERO)
   M%AVG_DROP_DEN=0._EB
   ALLOCATE(M%AVG_DROP_AREA(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_AREA',IZERO)
   M%AVG_DROP_AREA=0._EB
   ALLOCATE(M%AVG_DROP_TMP(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_TMP',IZERO)
   M%AVG_DROP_TMP=TMPM
   ALLOCATE(M%AVG_DROP_RAD(0:IBP1,0:JBP1,0:KBP1,N_LP_ARRAY_INDICES),STAT=IZERO)
   CALL ChkMemErr('INIT','AVG_DROP_RAD',IZERO)
   M%AVG_DROP_RAD=0._EB
ENDIF

IF (N_LP_ARRAY_INDICES>0 .OR. N_REACTIONS>0 .OR. ANY(SPECIES_MIXTURE%DEPOSITING) .OR. &
    ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0) .OR. REACTING_THIN_OBSTRUCTIONS .OR. INCLUDE_PYROLYSIS) THEN
   ALLOCATE(M%D_SOURCE(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','D_SOURCE',IZERO)
   M%D_SOURCE = 0._EB
   ALLOCATE(M%M_DOT_PPP(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES),STAT=IZERO)
   CALL ChkMemErr('INIT','M_DOT_PPP',IZERO)
   M%M_DOT_PPP=0._EB
ENDIF

! If radiation absorption desired allocate arrays

ALLOCATE(M%CHI_R(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)     ; CALL ChkMemErr('INIT','CHI_R',IZERO)     ; M%CHI_R     = 0._EB
ALLOCATE(M%QR(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)        ; CALL ChkMemErr('INIT','QR',IZERO)        ; M%QR        = 0._EB
ALLOCATE(M%KAPPA_GAS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','KAPPA_GAS',IZERO) ; M%KAPPA_GAS = 0._EB
ALLOCATE(M%UII(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)       ; CALL ChkMemErr('INIT','UII',IZERO)       ; M%UII       = 0._EB

! Work arrays

ALLOCATE(M%WORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK1',IZERO)
ALLOCATE(M%WORK2(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK2',IZERO)
ALLOCATE(M%WORK3(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK3',IZERO)
ALLOCATE(M%WORK4(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK4',IZERO)
ALLOCATE(M%WORK5(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK5',IZERO)
ALLOCATE(M%WORK6(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK6',IZERO)
ALLOCATE(M%WORK7(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK7',IZERO)
ALLOCATE(M%WORK8(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK8',IZERO)
ALLOCATE(M%WORK9(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','WORK9',IZERO)

ALLOCATE(M%IWORK1(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)                         ; CALL ChkMemErr('INIT','IWORK1',IZERO)
ALLOCATE(M%SCALAR_WORK1(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK1',IZERO)
ALLOCATE(M%SCALAR_WORK2(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK2',IZERO)
ALLOCATE(M%SCALAR_WORK3(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK3',IZERO)
ALLOCATE(M%SCALAR_WORK4(0:IBP1,0:JBP1,0:KBP1,1:N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','SCALAR_WORK4',IZERO)
M%IWORK1=0
M%SCALAR_WORK1=0._EB
M%SCALAR_WORK2=0._EB
M%SCALAR_WORK3=0._EB
M%SCALAR_WORK4=0._EB

IF (STRATIFICATION) THEN

   ! Compute the atmospheric pressure profile ramp using the specified temperature ramp

   RP => RAMPS(I_RAMP_P0_Z)
   INTEGRAL = 0._EB
   IF (HVAC_SOLVE) THEN
      ZSW = MIN(ZS_MIN-DZS_MAX,NODE_Z_MIN)
   ELSE
      ZSW = ZS_MIN
   ENDIF

   DO K=0,RP%NUMBER_INTERPOLATION_POINTS+1
      TEMP = TMPA*RAMPS(I_RAMP_TMP0_Z)%INTERPOLATED_DATA(K)
      INTEGRAL = INTEGRAL + (GVEC(3)/(RSUM0*TEMP))/RP%RDT
      RP%INTERPOLATED_DATA(K) = P_INF*EXP(GVEC(3)*(ZSW-GROUND_LEVEL)/(RSUM0*TMPA))*EXP(INTEGRAL)
   ENDDO

   ! Populate the cell-centered background temperature and pressure

   DO K=0,M%KBP1
      M%TMP_0(K) = TMPA*EVALUATE_RAMP(M%ZC(K),I_RAMP_TMP0_Z)
      M%P_0(K)   = EVALUATE_RAMP(M%ZC(K),I_RAMP_P0_Z)
   ENDDO
ELSE

   M%TMP_0(:) = TMPA
   M%P_0(:)   = P_INF

ENDIF

! Initialize density profile

DO K=0,M%KBP1
   M%RHO_0(K) = M%P_0(K)/(M%TMP_0(K)*RSUM0)
ENDDO

! Initialize various time step variables

DT_INITIAL = DT

! Initialize major arrays

ALLOCATE(M%U_WIND(0:M%KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','U_WIND',IZERO)
ALLOCATE(M%V_WIND(0:M%KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','V_WIND',IZERO)
ALLOCATE(M%W_WIND(0:M%KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','W_WIND',IZERO)

CALL COMPUTE_WIND_COMPONENTS(T_BEGIN,NM)

DO K=0,M%KBP1
   M%RHO(:,:,K)  = M%RHO_0(K)
   M%RHOS(:,:,K) = M%RHO_0(K)
   M%TMP(:,:,K)  = M%TMP_0(K)
   M%U(:,:,K)    = M%U_WIND(K)
   M%V(:,:,K)    = M%V_WIND(K)
   M%W(:,:,K)    = M%W_WIND(K)
ENDDO

! Custom velocity RAMPS (for verification)

IF (I_RAMP_UX>0) THEN
   DO I=0,M%IBAR
      M%U(I,:,:) = EVALUATE_RAMP(M%X(I),I_RAMP_UX)
   ENDDO
ENDIF
IF (I_RAMP_UY>0) THEN
   DO J=0,M%JBP1
      M%U(:,J,:) = EVALUATE_RAMP(M%YC(J),I_RAMP_UY)
   ENDDO
ENDIF
IF (I_RAMP_UZ>0) THEN
   DO K=0,M%KBP1
      M%U(:,:,K) = EVALUATE_RAMP(M%ZC(K),I_RAMP_UZ)
   ENDDO
ENDIF

IF (I_RAMP_VX>0) THEN
   DO I=0,M%IBP1
      M%V(I,:,:) = EVALUATE_RAMP(M%XC(I),I_RAMP_VX)
   ENDDO
ENDIF
IF (I_RAMP_VY>0) THEN
   DO J=0,M%JBAR
      M%V(:,J,:) = EVALUATE_RAMP(M%Y(J),I_RAMP_VY)
   ENDDO
ENDIF
IF (I_RAMP_VZ>0) THEN
   DO K=0,M%KBP1
      M%V(:,:,K) = EVALUATE_RAMP(M%ZC(K),I_RAMP_VZ)
   ENDDO
ENDIF

IF (I_RAMP_WX>0) THEN
   DO I=0,M%IBP1
      M%W(I,:,:) = EVALUATE_RAMP(M%XC(I),I_RAMP_WX)
   ENDDO
ENDIF
IF (I_RAMP_WY>0) THEN
   DO J=0,M%JBP1
      M%W(:,J,:) = EVALUATE_RAMP(M%YC(J),I_RAMP_WY)
   ENDDO
ENDIF
IF (I_RAMP_WZ>0) THEN
   DO K=0,M%KBAR
      M%W(:,:,K) = EVALUATE_RAMP(M%Z(K),I_RAMP_WZ)
   ENDDO
ENDIF

M%US    = M%U
M%VS    = M%V
M%WS    = M%W
M%FVX   = 0._EB
M%FVY   = 0._EB
M%FVZ   = 0._EB
M%KRES  = 0._EB
M%H     = 0.5_EB*(U0**2+V0**2+W0**2)
M%HS    = 0.5_EB*(U0**2+V0**2+W0**2)
M%DDDT  = 0._EB
M%D     = 0._EB
M%DS    = 0._EB
M%Q     = 0._EB

! Calculate LES filter width

ALLOCATE(M%LES_FILTER_WIDTH(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','LES_FILTER_WIDTH',IZERO)

DO K=0,KBP1
   DO J=0,JBP1
      DO I=0,IBP1
         M%LES_FILTER_WIDTH(I,J,K) = LES_FILTER_WIDTH_FUNCTION(M%DX(I),M%DY(J),M%DZ(K))
      ENDDO
   ENDDO
ENDDO

! Viscosity

ZZ_GET(1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
CALL GET_VISCOSITY(ZZ_GET,MU_N,TMPA)
M%MU = MU_N

CS = C_SMAGORINSKY
SELECT CASE(TURB_MODEL)
   CASE(CONSMAG,DYNSMAG)
      DO K=0,KBP1
         DO J=0,JBP1
            DO I=0,IBP1
               DELTA = M%LES_FILTER_WIDTH(I,J,K)
               M%CSD2(I,J,K) = (CS*DELTA)**2
            ENDDO
         ENDDO
      ENDDO
END SELECT

! Initialize mass fraction arrays

DO N=1,N_TRACKED_SPECIES
   M%ZZ(:,:,:,N)  = SPECIES_MIXTURE(N)%ZZ0
   M%ZZS(:,:,:,N) = SPECIES_MIXTURE(N)%ZZ0
ENDDO
DO N=N_TRACKED_SPECIES+1,N_TOTAL_SCALARS
   M%ZZ(:,:,:,N)  = INITIAL_UNMIXED_FRACTION
   M%ZZS(:,:,:,N) = INITIAL_UNMIXED_FRACTION
ENDDO

! Allocate and Initialize Mesh-Dependent Radiation Arrays

M%QR  = 0._EB
M%UII = 4._EB*SIGMA*TMPA4
M%ANGLE_INC_COUNTER = 0
M%RAD_CALL_COUNTER  = 0
IF (RADIATION) THEN
   ALLOCATE(M%UIID(0:M%IBP1,0:M%JBP1,0:M%KBP1,1:UIIDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','UIID',IZERO)
   M%UIID = 4._EB*SIGMA*TMPA4/REAL(UIIDIM,EB)
ENDIF

! Over-ride default ambient gas species mass fractions, temperature and density

DO N=1,N_INIT
   IN => INITIALIZATION(N)
   IF ((IN%NODE_ID/='null')) CYCLE
   IF (.NOT. (IN%ADJUST_INITIAL_CONDITIONS)) CYCLE
   DO K=0,KBP1
      DO J=0,JBP1
         DO I=0,IBP1
            IF (M%XC(I)<IN%X1.OR.M%XC(I)>IN%X2.OR.M%YC(J)<IN%Y1.OR.M%YC(J)>IN%Y2.OR.M%ZC(K)<IN%Z1.OR.M%ZC(K)>IN%Z2) CYCLE
            IF (IN%VOLUME_FRACTIONS_SPECIFIED) THEN
               VF = 0._EB
               DO NS=2,N_TRACKED_SPECIES
                  IF (IN%RAMP_VF_Z_INDEX(NS)>0) THEN
                     VF(NS) = EVALUATE_RAMP(M%ZC(K),IN%RAMP_VF_Z_INDEX(NS))
                  ELSE
                     VF(NS) = IN%VOLUME_FRACTION(NS)
                  ENDIF
               ENDDO
               VF(1) = 1._EB - SUM(VF)
               M%ZZ(I,J,K,1:N_TRACKED_SPECIES) = VF(1:N_TRACKED_SPECIES)*SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%MW / &
                                             SUM(VF(1:N_TRACKED_SPECIES)*SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%MW)
            ELSEIF (IN%MASS_FRACTIONS_SPECIFIED) THEN
               DO NS=2,N_TRACKED_SPECIES
                  IF (IN%RAMP_MF_Z_INDEX(NS)>0) THEN
                     M%ZZ(I,J,K,NS) = EVALUATE_RAMP(M%ZC(K),IN%RAMP_MF_Z_INDEX(NS))
                  ELSE
                     M%ZZ(I,J,K,NS) = IN%MASS_FRACTION(NS)
                  ENDIF
               ENDDO
               M%ZZ(I,J,K,1) = 1._EB - SUM(M%ZZ(I,J,K,2:N_TRACKED_SPECIES))
            ENDIF
            M%ZZS(I,J,K,1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES)
            IF (IN%RAMP_TMP_Z_INDEX>0) THEN
               M%TMP(I,J,K) = TMPM + EVALUATE_RAMP(M%ZC(K),IN%RAMP_TMP_Z_INDEX)
            ELSEIF (IN%TEMPERATURE>0._EB) THEN
               M%TMP(I,J,K) = IN%TEMPERATURE
            ENDIF
            ZZ_GET(1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,M%RSUM(I,J,K))
            M%RHO(I,J,K)  = M%P_0(K)/(M%TMP(I,J,K)*M%RSUM(I,J,K))
            M%RHOS(I,J,K) = M%RHO(I,J,K)
            IF (RADIATION) THEN
               M%UII(I,J,K) = 4._EB*SIGMA*M%TMP(I,J,K)**4
               M%UIID(I,J,K,1:UIIDIM) = M%UII(I,J,K)/REAL(UIIDIM,EB)
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

! General work arrays

M%WORK1 = 0._EB
M%WORK2 = 0._EB
M%WORK3 = 0._EB
M%WORK4 = 0._EB
M%WORK5 = 0._EB
M%WORK6 = 0._EB
M%WORK7 = 0._EB

! Allocate lagrangian particle storage array and compute the dimensions of its components

M%NLP = 0
M%NLPDIM = 50
M%PARTICLE_TAG = NM
IF (N_LAGRANGIAN_CLASSES > 0) THEN
   ALLOCATE(M%PARTICLE_LAST(N_LAGRANGIAN_CLASSES))
   M%PARTICLE_LAST = 0
ENDIF

IF (PARTICLE_FILE) THEN
   ALLOCATE(M%LAGRANGIAN_PARTICLE(M%NLPDIM),STAT=IZERO)
   CALL ChkMemErr('INIT','PARTICLE',IZERO)
ENDIF

! Allocate wall cell and the various BOUNDARY arrays. These arrays will grow as needed.

N_EXTERNAL_CELLS = 2*(M%IBAR*M%JBAR+M%IBAR*M%KBAR+M%JBAR*M%KBAR)

M%N_WALL_CELLS_DIM      = N_EXTERNAL_CELLS
M%N_THIN_WALL_CELLS_DIM = 10

ALLOCATE(M%WALL(0:M%N_WALL_CELLS_DIM),STAT=IZERO)           ; CALL ChkMemErr('INIT','WALL',IZERO)
ALLOCATE(M%THIN_WALL(0:M%N_THIN_WALL_CELLS_DIM),STAT=IZERO) ; CALL ChkMemErr('INIT','THIN_WALL',IZERO)

M%WALL(0)%BOUNDARY_TYPE = NULL_BOUNDARY
M%WALL(0)%SURF_INDEX    = DEFAULT_SURF_INDEX

! Allocate arrays that are dimensioned by the number of external wall cells

ALLOCATE(M%UVW_SAVE(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','UVW_SAVE',IZERO)
M%UVW_SAVE = 0._EB

ALLOCATE(M%U_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','U_GHOST',IZERO)
ALLOCATE(M%V_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','V_GHOST',IZERO)
ALLOCATE(M%W_GHOST(M%N_EXTERNAL_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','W_GHOST',IZERO)
M%U_GHOST = 0._EB
M%V_GHOST = 0._EB
M%W_GHOST = 0._EB

! Allocate arrays for turbulent inflow boundary conditions (experimental)

VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   EDDY_IF: IF (VT%N_EDDY>0) THEN
      SELECT CASE(ABS(VT%IOR))
         CASE(1)
            ALLOCATE(VT%U_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%J1+1:VT%J2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
         CASE(2)
            ALLOCATE(VT%U_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%I1+1:VT%I2,VT%K1+1:VT%K2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
         CASE(3)
            ALLOCATE(VT%U_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','U_EDDY',IZERO)
            ALLOCATE(VT%V_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','V_EDDY',IZERO)
            ALLOCATE(VT%W_EDDY(VT%I1+1:VT%I2,VT%J1+1:VT%J2),STAT=IZERO)
            CALL ChkMemErr('READ_VENT','W_EDDY',IZERO)
      END SELECT
      ALLOCATE(VT%X_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','X_EDDY',IZERO)
      ALLOCATE(VT%Y_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','Y_EDDY',IZERO)
      ALLOCATE(VT%Z_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','Z_EDDY',IZERO)
      ALLOCATE(VT%CU_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CU_EDDY',IZERO)
      ALLOCATE(VT%CV_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CV_EDDY',IZERO)
      ALLOCATE(VT%CW_EDDY(VT%N_EDDY),STAT=IZERO)
      CALL ChkMemErr('READ_VENT','CW_EDDY',IZERO)
      VT%U_EDDY=0._EB
      VT%V_EDDY=0._EB
      VT%W_EDDY=0._EB
      VT%X_EDDY=0._EB
      VT%Y_EDDY=0._EB
      VT%Z_EDDY=0._EB
      VT%CU_EDDY=0._EB
      VT%CV_EDDY=0._EB
      VT%CW_EDDY=0._EB
   ENDIF EDDY_IF
ENDDO VENT_LOOP

! Set up WALL for external boundaries of the current mesh

M%N_WALL_CELLS = 0

DO IOR=1,-1,-2
   IF (IOR== 1) I = 0
   IF (IOR==-1) I = IBP1
   DO K=1,KBAR
      DO J=1,JBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO

DO IOR=2,-2,-4
   IF (IOR== 2) J = 0
   IF (IOR==-2) J = JBP1
   DO K=1,KBAR
      DO I=1,IBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO

DO IOR=3,-3,-6
   IF (IOR== 3) K = 0
   IF (IOR==-3) K = KBP1
   DO J=1,JBAR
      DO I=1,IBAR
         M%N_WALL_CELLS = M%N_WALL_CELLS + 1
         CALL INIT_WALL_CELL(NM,I,J,K,0,M%N_WALL_CELLS,IOR,DEFAULT_SURF_INDEX,IERR,T_BEGIN) ; IF (IERR>0) RETURN
      ENDDO
   ENDDO
ENDDO

! Go through all obstructions and decide which cell faces ought to be given a wall cell index and initialized

M%N_INTERNAL_WALL_CELLS = 0
M%N_THIN_WALL_CELLS = 0

OBST_LOOP_2: DO OBST_INDEX=1,M%N_OBST

   OB=>M%OBSTRUCTION(OBST_INDEX)

   IF (ANY(SURFACE(OB%SURF_INDEX(:))%HT_DIM>1) .AND. (OB%I1==OB%I2 .OR. OB%J1==OB%J2 .OR. OB%K1==OB%K2)) THEN

      IF (OB%I1==OB%I2 .AND. ABS(OB%X2-OB%X1)>TWO_EPSILON_EB .AND. OB%UNDIVIDED_INPUT_LENGTH(1)<0.5_EB*M%DX(OB%I1)) THEN
         DO K=OB%K1+1,OB%K2
            IF (OB%J1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,OB%J1,K,OBST_INDEX,M%N_THIN_WALL_CELLS,-2,OB%SURF_INDEX(-2),3)
            ENDIF
            IF (OB%J2<M%JBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,OB%J2,K,OBST_INDEX,M%N_THIN_WALL_CELLS, 2,OB%SURF_INDEX( 2),3)
            ENDIF
         ENDDO
         DO J=OB%J1+1,OB%J2
            IF (OB%K1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,J,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS,-3,OB%SURF_INDEX(-3),2)
            ENDIF
            IF (OB%K2<M%KBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,J,OB%K2,OBST_INDEX,M%N_THIN_WALL_CELLS, 3,OB%SURF_INDEX( 3),2)
            ENDIF
         ENDDO
      ENDIF

      IF (OB%J1==OB%J2 .AND. ABS(OB%Y2-OB%Y1)>TWO_EPSILON_EB .AND. OB%UNDIVIDED_INPUT_LENGTH(2)<0.5_EB*M%DY(OB%J1)) THEN
         DO K=OB%K1+1,OB%K2
            IF (OB%I1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,OB%J1,K,OBST_INDEX,M%N_THIN_WALL_CELLS,-1,OB%SURF_INDEX(-1),3)
            ENDIF
            IF (OB%I2<M%IBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I2,OB%J1,K,OBST_INDEX,M%N_THIN_WALL_CELLS, 1,OB%SURF_INDEX( 1),3)
            ENDIF
         ENDDO
         DO I=OB%I1+1,OB%I2
            IF (OB%K1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,I,OB%J1,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS,-3,OB%SURF_INDEX(-3),1)
            ENDIF
            IF (OB%K2<M%KBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,I,OB%J1,OB%K2,OBST_INDEX,M%N_THIN_WALL_CELLS, 3,OB%SURF_INDEX( 3),1)
            ENDIF
         ENDDO
      ENDIF

      IF (OB%K1==OB%K2 .AND. ABS(OB%Z2-OB%Z1)>TWO_EPSILON_EB .AND. OB%UNDIVIDED_INPUT_LENGTH(3)<0.5_EB*M%DZ(OB%K1)) THEN
         DO I=OB%I1+1,OB%I2
            IF (OB%J1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,I,OB%J1,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS,-2,OB%SURF_INDEX(-2),1)
            ENDIF
            IF (OB%J2<M%JBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,I,OB%J2,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS, 2,OB%SURF_INDEX( 2),1)
            ENDIF
         ENDDO
         DO J=OB%J1+1,OB%J2
            IF (OB%I1>0) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I1,J,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS,-1,OB%SURF_INDEX(-1),2)
            ENDIF
            IF (OB%I2<M%IBAR) THEN
               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
               CALL INIT_THIN_WALL_CELL(NM,OB%I2,J,OB%K1,OBST_INDEX,M%N_THIN_WALL_CELLS, 1,OB%SURF_INDEX( 1),2)
            ENDIF
         ENDDO
      ENDIF

   ENDIF

   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (I==1) THEN ; OB%EXPOSED_FACE_INDEX(1)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I-1,J,K)
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE ! Permanently covered face
         IOR = -1
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO J=OB%J1+1,OB%J2
         I = OB%I2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (I==M%IBAR) THEN ; OB%EXPOSED_FACE_INDEX(2)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I+1,J,K)
         ! Permanently covered face
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE
         IOR = 1
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (J==1) THEN ; OB%EXPOSED_FACE_INDEX(3)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I,J-1,K)
         ! Permanently covered face
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE
         IOR = -2
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO K=OB%K1+1,OB%K2
      DO I=OB%I1+1,OB%I2
         J = OB%J2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (J==M%JBAR) THEN ; OB%EXPOSED_FACE_INDEX(4)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I,J+1,K)
         ! Permanently covered face
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE
         IOR = 2
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K1+1
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (K==1) THEN ; OB%EXPOSED_FACE_INDEX(5)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I,J,K-1)
         ! Permanently covered face
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE
         IOR = -3
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

   DO J=OB%J1+1,OB%J2
      DO I=OB%I1+1,OB%I2
         K = OB%K2
         ! Don't assign wall cell index to obstruction face pointing out of the computational domain
         IF (K==M%KBAR) THEN ; OB%EXPOSED_FACE_INDEX(6)=1 ; CYCLE ; ENDIF
         IC = M%CELL_INDEX(I,J,K+1)
         ! Permanently covered face
         IF (M%CELL(IC)%SOLID .AND. .NOT.M%OBSTRUCTION(M%CELL(IC)%OBST_INDEX)%REMOVABLE) CYCLE
         IOR = 3
         SURF_INDEX = OB%SURF_INDEX(IOR)
         IW  = M%CELL(IC)%WALL_INDEX(-IOR)
         IF (IW==0) THEN
            M%N_INTERNAL_WALL_CELLS = M%N_INTERNAL_WALL_CELLS + 1
            M%N_WALL_CELLS = M%N_EXTERNAL_WALL_CELLS + M%N_INTERNAL_WALL_CELLS
            IW  = M%N_WALL_CELLS
         ELSE
            IF (.NOT.OB%OVERLAY .OR. OB%HIDDEN) CYCLE
            CALL CHECK_OVERLAPPING_OBSTRUCTIONS(IERR) ; IF (IERR>0) RETURN
         ENDIF
         CALL INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,T_BEGIN)
         IF (IERR>0) RETURN
      ENDDO
   ENDDO

ENDDO OBST_LOOP_2

! For MULT/OBST/SHAPE, recompute B1%AREA_ADJUST

OBST_SHAPE_IF: IF (OBST_SHAPE_AREA_ADJUST) THEN

   ! First, sum the face areas of the OBSTs with a given SURF
   SHAPE_LOOP_1: DO IW=M%N_EXTERNAL_WALL_CELLS+1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
      WC=>M%WALL(IW)
      IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE SHAPE_LOOP_1
      BC=>M%BOUNDARY_COORD(WC%BC_INDEX)
      B1=>M%BOUNDARY_PROP1(WC%B1_INDEX)
      II=BC%II
      JJ=BC%JJ
      KK=BC%KK
      IOR=BC%IOR
      IC=M%CELL_INDEX(II,JJ,KK)
      OBST_INDEX=M%CELL(IC)%OBST_INDEX; IF (OBST_INDEX==0)       CYCLE SHAPE_LOOP_1
      OB=>M%OBSTRUCTION(OBST_INDEX); IF (OB%MULT_INDEX<0) CYCLE SHAPE_LOOP_1
      MR=>MULTIPLIER(OB%MULT_INDEX)
      SHAPE_SELECT_1: SELECT CASE(OB%SHAPE_TYPE)
         CASE(OBST_SPHERE_TYPE)
            MR%FDS_AREA(1) = MR%FDS_AREA(1) + B1%AREA
         CASE(OBST_CYLINDER_TYPE)
            ! OB%SHAPE_AREA follows the same pattern as SURF_IDS: top, sides, bottom
            SELECT CASE(IOR)
               CASE(3);         MR%FDS_AREA(1) = MR%FDS_AREA(1) + B1%AREA ! top
               CASE(-1,1,-2,2); MR%FDS_AREA(2) = MR%FDS_AREA(2) + B1%AREA ! side
               CASE(-3);        MR%FDS_AREA(3) = MR%FDS_AREA(3) + B1%AREA ! bottom
            END SELECT
         CASE(OBST_CONE_TYPE)
            SELECT CASE(IOR)
               CASE(-1,1,-2,2,3); MR%FDS_AREA(1) = MR%FDS_AREA(1) + B1%AREA
               CASE(-3);          MR%FDS_AREA(2) = MR%FDS_AREA(2) + B1%AREA
            END SELECT
         CASE(OBST_BOX_TYPE)
            ! Follows sextuplet ordering from SURF_ID6
            SELECT CASE(IOR)
               CASE(-1); MR%FDS_AREA(1) = MR%FDS_AREA(1) + B1%AREA
               CASE( 1); MR%FDS_AREA(2) = MR%FDS_AREA(2) + B1%AREA
               CASE(-2); MR%FDS_AREA(3) = MR%FDS_AREA(3) + B1%AREA
               CASE( 2); MR%FDS_AREA(4) = MR%FDS_AREA(4) + B1%AREA
               CASE(-3); MR%FDS_AREA(5) = MR%FDS_AREA(5) + B1%AREA
               CASE( 3); MR%FDS_AREA(6) = MR%FDS_AREA(6) + B1%AREA
            END SELECT
      END SELECT SHAPE_SELECT_1
   ENDDO SHAPE_LOOP_1

   ! Next, AREA_ADJUST the WALL_CELL

   SHAPE_LOOP_2: DO IW=M%N_EXTERNAL_WALL_CELLS+1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
      WC=>M%WALL(IW)
      IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE SHAPE_LOOP_2
      BC=>M%BOUNDARY_COORD(WC%BC_INDEX)
      B1=>M%BOUNDARY_PROP1(WC%B1_INDEX)
      II=BC%II
      JJ=BC%JJ
      KK=BC%KK
      IOR=BC%IOR
      IC=M%CELL_INDEX(II,JJ,KK)
      OBST_INDEX=M%CELL(IC)%OBST_INDEX; IF (OBST_INDEX==0)       CYCLE SHAPE_LOOP_2
      OB=>M%OBSTRUCTION(OBST_INDEX); IF (OB%MULT_INDEX<0) CYCLE SHAPE_LOOP_2
      MR=>MULTIPLIER(OB%MULT_INDEX)
      SF=>SURFACE(WC%SURF_INDEX)
      IF (OB%SHAPE_TYPE>0) B1%AREA_ADJUST=1._EB
      SHAPE_SELECT_2: SELECT CASE(OB%SHAPE_TYPE)
         CASE(OBST_SPHERE_TYPE)
            B1%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
         CASE(OBST_CYLINDER_TYPE)
            SELECT CASE(IOR)
               CASE(3);         B1%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1) ! top
               CASE(-1,1,-2,2); B1%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2) ! side
               CASE(-3);        B1%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(3) ! bottom
            END SELECT
         CASE(OBST_CONE_TYPE)
            SELECT CASE(IOR)
               CASE(-1,1,-2,2,3); B1%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
               CASE(-3);          B1%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2)
            END SELECT
         CASE(OBST_BOX_TYPE)
            SELECT CASE(IOR)
               CASE(-1); B1%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(1)
               CASE( 1); B1%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(2)
               CASE(-2); B1%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(3)
               CASE( 2); B1%AREA_ADJUST = OB%SHAPE_AREA(1)/MR%FDS_AREA(4)
               CASE(-3); B1%AREA_ADJUST = OB%SHAPE_AREA(2)/MR%FDS_AREA(5)
               CASE( 3); B1%AREA_ADJUST = OB%SHAPE_AREA(3)/MR%FDS_AREA(6)
            END SELECT
      END SELECT SHAPE_SELECT_2
      B1%AREA_ADJUST = B1%AREA_ADJUST*SF%AREA_MULTIPLIER
   ENDDO SHAPE_LOOP_2

ENDIF OBST_SHAPE_IF

! Reset ghost cell values of cell centered velocity for use in computing viscosity (must be done after INIT_WALL_CELL)

DO IW=1,M%N_EXTERNAL_WALL_CELLS
   WC=>M%WALL(IW)
   BC=>M%BOUNDARY_COORD(WC%BC_INDEX)
   M%U_GHOST(IW) = M%U_WIND(BC%KKG)
   M%V_GHOST(IW) = M%V_WIND(BC%KKG)
   M%W_GHOST(IW) = M%W_WIND(BC%KKG)
ENDDO

CONTAINS


!> \brief Check if two removable obstructions overlap at a common surface.
!> \param IERR Error flag.

SUBROUTINE CHECK_OVERLAPPING_OBSTRUCTIONS(IERR)

USE COMP_FUNCTIONS, ONLY: SHUTDOWN
INTEGER, INTENT(OUT) :: IERR
INTEGER :: OBST_INDEX_PREVIOUS,OB_SURF_INDEX,OB_PREVIOUS_SURF_INDEX
TYPE(OBSTRUCTION_TYPE), POINTER :: OB_PREVIOUS

IERR = 0
RETURN
OBST_INDEX_PREVIOUS = M%WALL(IW)%OBST_INDEX
OB_PREVIOUS => M%OBSTRUCTION(OBST_INDEX_PREVIOUS)
OB_SURF_INDEX = OB%SURF_INDEX(IOR)
OB_PREVIOUS_SURF_INDEX = OB_PREVIOUS%SURF_INDEX(IOR)
 
IF (OBST_INDEX_PREVIOUS>0 .AND. OBST_INDEX_PREVIOUS/=OBST_INDEX) THEN
   IF ( (OB%REMOVABLE.OR.OB_PREVIOUS%REMOVABLE) .AND. OB_SURF_INDEX/=OB_PREVIOUS_SURF_INDEX .AND. &
        (SURFACE(OB_SURF_INDEX)%THERMAL_BC_INDEX==THERMALLY_THICK .OR.  &
         SURFACE(OB_PREVIOUS_SURF_INDEX)%THERMAL_BC_INDEX==THERMALLY_THICK)) THEN
      WRITE(LU_ERR,'(5A,I0)') 'WARNING(613): OBST ',TRIM(OB%ID),' and OBST ',TRIM(OB_PREVIOUS%ID),' overlap surfaces in Mesh ',NM
      IERR = 0
   ENDIF
ENDIF

END SUBROUTINE CHECK_OVERLAPPING_OBSTRUCTIONS

END SUBROUTINE INITIALIZE_MESH_VARIABLES_1


!> \brief Continuation of variable allocation and other setup functions
!> \param NM Mesh number

SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM)

USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_EDGE,REALLOCATE_REAL_ARRAY
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_HEAT
USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE CONTROL_VARIABLES
INTEGER :: N,I,J,K,IPTS,JPTS,KPTS,N_EDGES_DIM,IW,IC,IERR,IPZ,IZERO,ICF,NSLICE
INTEGER, INTENT(IN) :: NM
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),VC,RTRM,CP
INTEGER :: IBP1,JBP1,KBP1,IBAR,JBAR,KBAR
REAL(EB) :: XS,XF,YS,YF,ZS,ZF
TYPE (MESH_TYPE), POINTER :: M
TYPE (WALL_TYPE), POINTER :: WC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE (CFACE_TYPE), POINTER :: CFA
TYPE (MESH_TYPE), POINTER :: OM
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (SURFACE_TYPE), POINTER :: SF
LOGICAL :: SOLID_CELL

IERR = 0
M => MESHES(NM)
IBP1 = M%IBP1
JBP1 = M%JBP1
KBP1 = M%KBP1
IBAR = M%IBAR
JBAR = M%JBAR
KBAR = M%KBAR
XS = M%XS
YS = M%YS
ZS = M%ZS
XF = M%XF
YF = M%YF
ZF = M%ZF

! Surface work arrays

ALLOCATE(M%WALL_WORK1(M%N_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK1',IZERO)
ALLOCATE(M%WALL_WORK2(M%N_WALL_CELLS),STAT=IZERO)
CALL ChkMemErr('INIT','WALL_WORK2',IZERO)

! Background pressure variables

ALLOCATE(  M%PBAR(0:KBP1,0:N_ZONE),STAT=IZERO)   ; CALL ChkMemErr('INIT','PBAR',IZERO)
ALLOCATE(  M%PBAR_S(0:KBP1,0:N_ZONE),STAT=IZERO) ; CALL ChkMemErr('INIT','PBAR_S',IZERO)
ALLOCATE(  M%R_PBAR(0:KBP1,0:N_ZONE),STAT=IZERO) ; CALL ChkMemErr('INIT','R_PBAR',IZERO)
ALLOCATE(  M%D_PBAR_DT(0:N_ZONE),STAT=IZERO)     ; CALL ChkMemErr('INIT','D_PBAR_DT',IZERO)    ; M%D_PBAR_DT = 0._EB
ALLOCATE(  M%D_PBAR_DT_S(0:N_ZONE),STAT=IZERO)   ; CALL ChkMemErr('INIT','D_PBAR_DT_S',IZERO)  ; M%D_PBAR_DT_S = 0._EB
ALLOCATE(  M%U_LEAK(0:N_ZONE),STAT=IZERO)        ; CALL ChkMemErr('INIT','U_LEAK',IZERO)       ; M%U_LEAK = 0._EB

DO K=0,M%KBP1
   M%PBAR(K,:)   = M%P_0(K)
   M%PBAR_S(K,:) = M%P_0(K)
ENDDO

! Initialize PSUM for zone cases

IF (N_ZONE > 0) THEN
   ZONE_LOOP: DO IPZ = 1,N_ZONE
      PSUM(IPZ) = 0._EB
      DO K=1,M%KBAR
         DO J=1,M%JBAR
            DO I=1,M%IBAR
               IF (M%PRESSURE_ZONE(I,J,K) /= IPZ) CYCLE
               IF (M%CELL(M%CELL_INDEX(I,J,K))%SOLID)  CYCLE
               VC   = M%DX(I)*M%RC(I)*M%DY(J)*M%DZ(K)
               ZZ_GET(1:N_TRACKED_SPECIES) = M%ZZ(I,J,K,1:N_TRACKED_SPECIES)
               CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,M%TMP(I,J,K))
               RTRM = M%RSUM(I,J,K)/(CP*M%PBAR(K,IPZ))
               PSUM(IPZ) = PSUM(IPZ) + VC*(1._EB/M%PBAR(K,IPZ)-RTRM)
            ENDDO
         ENDDO
      ENDDO
   ENDDO ZONE_LOOP
ENDIF

! Loop through WALL and CFACE cells and assign PRESSURE_ZONE. Also, check for
! inappropriate boundaries, such as thin obstructions that burn or blow.

WALL_LOOP_0: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS

   WC => M%WALL(IW)
   BC => M%BOUNDARY_COORD(WC%BC_INDEX)
   B1 => M%BOUNDARY_PROP1(WC%B1_INDEX)
   SF => SURFACE(WC%SURF_INDEX)

   B1%PRESSURE_ZONE = M%PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG)

   IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN
      IF (M%EXTERNAL_WALL(IW)%NOM>0) THEN
         OM => MESHES(M%EXTERNAL_WALL(IW)%NOM)
         IC = OM%CELL_INDEX(M%EXTERNAL_WALL(IW)%IIO_MIN,M%EXTERNAL_WALL(IW)%JJO_MIN,M%EXTERNAL_WALL(IW)%KKO_MIN)
         SOLID_CELL = OM%CELL(IC)%SOLID
      ELSE
         IC =  M%CELL_INDEX(BC%II,BC%JJ,BC%KK)
         SOLID_CELL = M%CELL(IC)%SOLID
      ENDIF
   ELSE
      IC =  M%CELL_INDEX(BC%II,BC%JJ,BC%KK)
      SOLID_CELL = M%CELL(IC)%SOLID
   ENDIF

   IF (.NOT.SOLID_CELL) THEN
      IF ( (ABS(B1%U_NORMAL_0)>TWO_EPSILON_EB .OR. ANY(SF%LEAK_PATH>=0)) .AND. WC%OBST_INDEX>0 ) THEN
         WRITE(LU_ERR,'(A,A,A,I0)') 'ERROR(421): SURF ',TRIM(SF%ID),' cannot be applied to a thin obstruction, OBST #',&
                                    M%OBSTRUCTION(WC%OBST_INDEX)%ORDINAL
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
      IF (WC%VENT_INDEX>0 .AND. WC%OBST_INDEX>0) THEN
         VT => M%VENTS(WC%VENT_INDEX)
         IF (VT%BOUNDARY_TYPE==HVAC_BOUNDARY) THEN
            WRITE(LU_ERR,'(A,A,A,I0)') 'ERROR(422): VENT ',TRIM(VT%ID),' cannot be applied to a thin obstruction, OBST #',&
                                    M%OBSTRUCTION(WC%OBST_INDEX)%ORDINAL
            STOP_STATUS = SETUP_STOP
            RETURN
         ENDIF
      ENDIF
   ENDIF

ENDDO WALL_LOOP_0

CFACE_LOOP_0: DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTWALL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
   CFA => M%CFACE(ICF)
   BC  => M%BOUNDARY_COORD(CFA%BC_INDEX)
   B1  => M%BOUNDARY_PROP1(CFA%B1_INDEX)
   B1%PRESSURE_ZONE = M%PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG)
ENDDO CFACE_LOOP_0

! If there is complex terrain using GEOM and Above Ground Level (AGL) slices,
! determine K index of gas phase quantities.

DO NSLICE = 1, M%N_TERRAIN_SLCF
   IF (CC_IBM) THEN
      DO ICF=1,M%N_CUTFACE_MESH
         IF (M%CUT_FACE(ICF)%STATUS/=2 .OR. M%CUT_FACE(ICF)%NFACE<1) CYCLE
         IW  = MAXLOC(M%CUT_FACE(ICF)%AREA(1:M%CUT_FACE(ICF)%NFACE),DIM=1)
         CFA => M%CFACE( M%CUT_FACE(ICF)%CFACE_INDEX(IW) )
         BC  => M%BOUNDARY_COORD(CFA%BC_INDEX)
         IF (BC%NVEC(KAXIS)>-TWO_EPSILON_EB .AND.  CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
            IF (BC%KKG > M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)) THEN
               M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE) = MIN( M%KBAR , M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)+BC%KKG )
            ENDIF
         ENDIF
      ENDDO
   ELSE
      DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
         WC => M%WALL(IW)
         BC => M%BOUNDARY_COORD(WC%BC_INDEX)
         IF (BC%IOR/=3 .OR. WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE
         M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE) = MIN( M%KBAR , M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)+BC%KKG )
      ENDDO
   ENDIF
ENDDO

! Set clocks and counters related to frequency of solid phase conduction updates

M%BC_CLOCK     = T_BEGIN

! Allocate arrays for storing velocity boundary condition info

N_EDGES_DIM = 4*(IBP1*JBP1+IBP1*KBP1+JBP1*KBP1)
DO N=1,M%N_OBST
   OB=>M%OBSTRUCTION(N)
   IPTS = OB%I2-OB%I1
   JPTS = OB%J2-OB%J1
   KPTS = OB%K2-OB%K1
   N_EDGES_DIM = N_EDGES_DIM + 4*(IPTS*JPTS+IPTS*KPTS+JPTS*KPTS)
ENDDO

CALL REALLOCATE_EDGE(NM,N_EDGES_DIM,N_EDGES_DIM)

! Allocate array to hold character strings for Smokeview file

M%N_STRINGS     =   0
M%N_STRINGS_MAX = 100
ALLOCATE(M%STRING(M%N_STRINGS_MAX),STAT=IZERO)
CALL ChkMemErr('INIT','STRING',IZERO)

! Set up arrays to hold velocity boundary condition info

CALL INITIALIZE_EDGES

! Initialize Pressure solver

IF (.NOT.FREEZE_VELOCITY) CALL INITIALIZE_POISSON_SOLVER(NM)

IF (IERR/=0) RETURN

! Initialize Mesh Exchange

CALL INITIALIZE_INTERPOLATION


CONTAINS


!> \brief Set up edge arrays for velocity boundary conditions

SUBROUTINE INITIALIZE_EDGES

EDGE_COUNT(NM) = 0

! Arguments for DEFINE_EDGE(I,J,K,IOR,IEC,NM,OBST_INDEX)

DO K=0,M%KBAR
   DO J=0,M%JBAR
      IF (J>0) CALL DEFINE_EDGE(     0,J,K, 1,2,NM,0)
      IF (J>0) CALL DEFINE_EDGE(M%IBAR,J,K,-1,2,NM,0)
      IF (K>0) CALL DEFINE_EDGE(     0,J,K, 1,3,NM,0)
      IF (K>0) CALL DEFINE_EDGE(M%IBAR,J,K,-1,3,NM,0)
   ENDDO
ENDDO
DO K=0,M%KBAR
   DO I=0,M%IBAR
      IF (I>0) CALL DEFINE_EDGE(I,     0,K, 2,1,NM,0)
      IF (I>0) CALL DEFINE_EDGE(I,M%JBAR,K,-2,1,NM,0)
      IF (K>0) CALL DEFINE_EDGE(I,     0,K, 2,3,NM,0)
      IF (K>0) CALL DEFINE_EDGE(I,M%JBAR,K,-2,3,NM,0)
   ENDDO
ENDDO
DO J=0,M%JBAR
   DO I=0,M%IBAR
      IF (I>0) CALL DEFINE_EDGE(I,J,     0, 3,1,NM,0)
      IF (I>0) CALL DEFINE_EDGE(I,J,M%KBAR,-3,1,NM,0)
      IF (J>0) CALL DEFINE_EDGE(I,J,     0, 3,2,NM,0)
      IF (J>0) CALL DEFINE_EDGE(I,J,M%KBAR,-3,2,NM,0)
   ENDDO
ENDDO

IF (IERR/=0) RETURN

OBST_LOOP_3: DO N=1,M%N_OBST
   OB => M%OBSTRUCTION(N)
   DO K=OB%K1,OB%K2
      DO J=OB%J1,OB%J2
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I1,J,K,-1,2,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I2,J,K, 1,2,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I1,J,K,-1,3,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I2,J,K, 1,3,NM,N)
      ENDDO
   ENDDO
   DO K=OB%K1,OB%K2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J1,K,-2,1,NM,N)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J2,K, 2,1,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J1,K,-2,3,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J2,K, 2,3,NM,N)
      ENDDO
   ENDDO
   DO J=OB%J1,OB%J2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K1,-3,1,NM,N)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K2, 3,1,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K1,-3,2,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K2, 3,2,NM,N)
      ENDDO
   ENDDO
ENDDO OBST_LOOP_3

END SUBROUTINE INITIALIZE_EDGES


!> \brief Assign parameters for a given cell edge
!> \param II Index of edge in the x direction
!> \param JJ Index of edge in the y direction
!> \param KK Index of edge in the z direction
!> \param IOR Orientation index of adjacent wall cell
!> \param IEC Orientation of edge; 1=x direction; 2=y; 3=z
!> \param NM Mesh number
!> \param OBST_INDEX Obstruction index of edge

SUBROUTINE DEFINE_EDGE(II,JJ,KK,IOR,IEC,NM,OBST_INDEX)

INTEGER, INTENT(IN) :: II,JJ,KK,IOR,IEC,NM
INTEGER :: NOM,ICMM,ICMP,ICPM,ICPP,OBST_INDEX,IE,IW,IIO,JJO,KKO,IW1,IW2
REAL(EB) :: XI,YJ,ZK
TYPE (MESH_TYPE), POINTER :: MM
TYPE (EDGE_TYPE), POINTER :: ED

IF (OBST_INDEX>0) OB=>M%OBSTRUCTION(OBST_INDEX)

! Find the wall cells on each side of the edge

IW1 = -1
IW2 = -1

EDGE_DIRECTION_1: SELECT CASE(IEC)
   CASE(1) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-2)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ,KK)  )%WALL_INDEX(2)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(2)
         CASE( 2)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ+1,KK)  )%WALL_INDEX(-2)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ+1,KK+1))%WALL_INDEX(-2)
         CASE(-3)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ  ,KK))%WALL_INDEX(3)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(3)
         CASE( 3)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ  ,KK+1))%WALL_INDEX(-3)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ+1,KK+1))%WALL_INDEX(-3)
      END SELECT
   CASE(2) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ,KK)  )%WALL_INDEX(1)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ,KK+1))%WALL_INDEX(1)
         CASE( 1)
            IW1 = M%CELL(M%CELL_INDEX(II+1,JJ,KK)  )%WALL_INDEX(-1)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ,KK+1))%WALL_INDEX(-1)
         CASE(-3)
            IW1 = M%CELL(M%CELL_INDEX(II  ,JJ,KK))%WALL_INDEX(3)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(3)
         CASE( 3)
            IW1 = M%CELL(M%CELL_INDEX(II  ,JJ,KK+1))%WALL_INDEX(-3)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ,KK+1))%WALL_INDEX(-3)
      END SELECT
   CASE(3) EDGE_DIRECTION_1
      SELECT CASE(IOR)
         CASE(-1)
            IW1 = M%CELL(M%CELL_INDEX(II,JJ  ,KK))%WALL_INDEX(1)
            IW2 = M%CELL(M%CELL_INDEX(II,JJ+1,KK))%WALL_INDEX(1)
         CASE( 1)
            IW1 = M%CELL(M%CELL_INDEX(II+1,JJ  ,KK))%WALL_INDEX(-1)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ+1,KK))%WALL_INDEX(-1)
         CASE(-2)
            IW1 = M%CELL(M%CELL_INDEX(II  ,JJ,KK))%WALL_INDEX(2)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ,KK))%WALL_INDEX(2)
         CASE( 2)
            IW1 = M%CELL(M%CELL_INDEX(II  ,JJ+1,KK))%WALL_INDEX(-2)
            IW2 = M%CELL(M%CELL_INDEX(II+1,JJ+1,KK))%WALL_INDEX(-2)
      END SELECT
END SELECT EDGE_DIRECTION_1

! Decide what to do based on whether or not adjacent tiles exist

IF (IW1==0 .AND. IW2==0) RETURN
IF (IW1> 0 .AND. IW2==0) IW = IW1
IF (IW1==0 .AND. IW2> 0) IW = IW2
IF (IW1> 0 .AND. IW2> 0) THEN
   IW = IW2
   IF (IW1<=M%N_EXTERNAL_WALL_CELLS) THEN
      IF (M%EXTERNAL_WALL(IW1)%NOM>0) IW = IW1
   ENDIF
   IF (IW2<=M%N_EXTERNAL_WALL_CELLS) THEN
      IF (M%EXTERNAL_WALL(IW2)%NOM>0) IW = IW2
   ENDIF
ENDIF

! Assign the Index of the Edge (IE) and add to the list

ICMM = M%CELL_INDEX(II,JJ,KK)
SELECT CASE(IEC)
   CASE(1)
      IE = M%CELL(ICMM)%EDGE_INDEX( 4)
   CASE(2)
      IE = M%CELL(ICMM)%EDGE_INDEX( 8)
   CASE(3)
      IE = M%CELL(ICMM)%EDGE_INDEX(12)
END SELECT

IF (IE==0) THEN
   EDGE_COUNT(NM) = EDGE_COUNT(NM) + 1
   IE = EDGE_COUNT(NM)
ENDIF

ED => M%EDGE(IE)

! Determine the wall index of the adjacent wall tile

NOM = 0
IIO = 0
JJO = 0
KKO = 0

IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN
   IF (M%EXTERNAL_WALL(IW)%NOM>0) THEN
      NOM = M%EXTERNAL_WALL(IW)%NOM
      IIO = M%EXTERNAL_WALL(IW)%IIO_MIN
      JJO = M%EXTERNAL_WALL(IW)%JJO_MIN
      KKO = M%EXTERNAL_WALL(IW)%KKO_MIN
   ENDIF
ENDIF

! Identify EDGEs that lie at the external edge of the mesh

IF ( (II==0    .AND. KK==0   ) .OR. &
     (II==0    .AND. KK==KBAR) .OR. &
     (II==IBAR .AND. KK==0   ) .OR. &
     (II==IBAR .AND. KK==KBAR) .OR. &
     (II==0    .AND. JJ==0   ) .OR. &
     (II==0    .AND. JJ==JBAR) .OR. &
     (II==IBAR .AND. JJ==0   ) .OR. &
     (II==IBAR .AND. JJ==JBAR) .OR. &
     (JJ==0    .AND. KK==0   ) .OR. &
     (JJ==0    .AND. KK==KBAR) .OR. &
     (JJ==JBAR .AND. KK==0   ) .OR. &
     (JJ==JBAR .AND. KK==KBAR) ) ED%EXTERNAL=.TRUE.

! Fill up EDGE array

ED%I = II
ED%J = JJ
ED%K = KK
ED%AXIS = IEC

EDGE_DIRECTION_2: SELECT CASE(IEC)

   CASE (1) EDGE_DIRECTION_2

      ICPM = M%CELL_INDEX(II,JJ+1,KK)
      ICPP = M%CELL_INDEX(II,JJ+1,KK+1)
      ICMP = M%CELL_INDEX(II,JJ,KK+1)
      ED%CELL_INDEX_MM = ICMM
      ED%CELL_INDEX_PM = ICPM
      ED%CELL_INDEX_MP = ICMP
      ED%CELL_INDEX_PP = ICPP
      M%CELL(ICPP)%EDGE_INDEX(1) = IE
      M%CELL(ICMP)%EDGE_INDEX(2) = IE
      M%CELL(ICPM)%EDGE_INDEX(3) = IE
      M%CELL(ICMM)%EDGE_INDEX(4) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE(2)
               IF (IOR>0) ED%NOM_1 = -NOM
               IF (IOR<0) ED%NOM_1 =  NOM
               ED%IIO_1 = IIO
               ED%JJO_1 = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ONE_M_EPS , MM%CELLSK(NINT((M%Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               ED%EDGE_INTERPOLATION_FACTOR(1) = ZK-KKO
               ED%KKO_1 = KKO
            CASE(3)
               IF (IOR>0) ED%NOM_2 = -NOM
               IF (IOR<0) ED%NOM_2 =  NOM
               ED%IIO_2 = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ONE_M_EPS , MM%CELLSJ(NINT((M%Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               ED%EDGE_INTERPOLATION_FACTOR(2) = YJ-JJO
               ED%JJO_2 = JJO
               ED%KKO_2 = KKO
         END SELECT
      ENDIF

   CASE (2) EDGE_DIRECTION_2

      ICMP = M%CELL_INDEX(II+1,JJ,KK)
      ICPP = M%CELL_INDEX(II+1,JJ,KK+1)
      ICPM = M%CELL_INDEX(II,JJ,KK+1)
      ED%CELL_INDEX_MM = ICMM
      ED%CELL_INDEX_PM = ICPM
      ED%CELL_INDEX_MP = ICMP
      ED%CELL_INDEX_PP = ICPP
      M%CELL(ICPP)%EDGE_INDEX(5) = IE
      M%CELL(ICPM)%EDGE_INDEX(6) = IE
      M%CELL(ICMP)%EDGE_INDEX(7) = IE
      M%CELL(ICMM)%EDGE_INDEX(8) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) ED%NOM_2 = -NOM
               IF (IOR<0) ED%NOM_2 =  NOM
               ED%IIO_2 = IIO
               ED%JJO_2 = JJO
               MM => MESHES(NOM)
               ZK  = MIN( REAL(MM%KBAR,EB)+ONE_M_EPS , MM%CELLSK(NINT((M%Z(KK)-MM%ZS)*MM%RDZINT))+1._EB )
               KKO = MAX(1,FLOOR(ZK))
               ED%EDGE_INTERPOLATION_FACTOR(2) = ZK-KKO
               ED%KKO_2 = KKO
            CASE( 3)
               IF (IOR>0) ED%NOM_1 = -NOM
               IF (IOR<0) ED%NOM_1 =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ONE_M_EPS , MM%CELLSI(NINT((M%X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               ED%EDGE_INTERPOLATION_FACTOR(1) = XI-IIO
               ED%IIO_1 = IIO
               ED%JJO_1 = JJO
               ED%KKO_1 = KKO
         END SELECT
      ENDIF

   CASE (3) EDGE_DIRECTION_2

      ICPM = M%CELL_INDEX(II+1,JJ,KK)
      ICPP = M%CELL_INDEX(II+1,JJ+1,KK)
      ICMP = M%CELL_INDEX(II,JJ+1,KK)
      ED%CELL_INDEX_MM = ICMM
      ED%CELL_INDEX_PM = ICPM
      ED%CELL_INDEX_MP = ICMP
      ED%CELL_INDEX_PP = ICPP
      M%CELL(ICPP)%EDGE_INDEX( 9) = IE
      M%CELL(ICMP)%EDGE_INDEX(10) = IE
      M%CELL(ICPM)%EDGE_INDEX(11) = IE
      M%CELL(ICMM)%EDGE_INDEX(12) = IE
      IF (NOM/=0) THEN
         SELECT CASE(ABS(IOR))
            CASE( 1)
               IF (IOR>0) ED%NOM_1 = -NOM
               IF (IOR<0) ED%NOM_1 =  NOM
               ED%IIO_1 = IIO
               MM => MESHES(NOM)
               YJ  = MIN( REAL(MM%JBAR,EB)+ONE_M_EPS , MM%CELLSJ(NINT((M%Y(JJ)-MM%YS)*MM%RDYINT))+1._EB )
               JJO = MAX(1,FLOOR(YJ))
               ED%EDGE_INTERPOLATION_FACTOR(1) = YJ-JJO
               ED%JJO_1 = JJO
               ED%KKO_1 = KKO
            CASE( 2)
               IF (IOR>0) ED%NOM_2 = -NOM
               IF (IOR<0) ED%NOM_2 =  NOM
               MM => MESHES(NOM)
               XI  = MIN( REAL(MM%IBAR,EB)+ONE_M_EPS , MM%CELLSI(NINT((M%X(II)-MM%XS)*MM%RDXINT))+1._EB )
               IIO = MAX(1,FLOOR(XI))
               ED%EDGE_INTERPOLATION_FACTOR(2) = XI-IIO
               ED%IIO_2 = IIO
               ED%JJO_2 = JJO
               ED%KKO_2 = KKO
         END SELECT
      ENDIF

END SELECT EDGE_DIRECTION_2

END SUBROUTINE DEFINE_EDGE


!> \brief Create arrays by which info is to exchanged across meshes

SUBROUTINE INITIALIZE_INTERPOLATION

INTEGER :: NOM,I,J,K
TYPE (MESH_TYPE), POINTER :: M2

ALLOCATE(M%INTERPOLATED_MESH(1:M%IBAR,1:M%JBAR,1:M%KBAR),  STAT=IZERO)
CALL ChkMemErr('INIT','INTERPOLATED_MESH',IZERO)
M%INTERPOLATED_MESH = 0

DO K=1,M%KBAR
   DO J=1,M%JBAR
      DO I=1,M%IBAR
         OTHER_MESH_LOOP: DO NOM=1,NM-1
            M2=>MESHES(NOM)
            IF (M%X(I-1)>=M2%XS .AND. M%X(I)<=M2%XF .AND. &
                M%Y(J-1)>=M2%YS .AND. M%Y(J)<=M2%YF .AND. &
                M%Z(K-1)>=M2%ZS .AND. M%Z(K)<=M2%ZF) THEN
               M%INTERPOLATED_MESH(I,J,K) = NOM
               EXIT OTHER_MESH_LOOP
            ENDIF
         ENDDO OTHER_MESH_LOOP
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE INITIALIZE_INTERPOLATION

END SUBROUTINE INITIALIZE_MESH_VARIABLES_2


!> \brief Find WALL THIN_WALL cells with VARIABLE_THICKNESS or HT3D and adjust the 1-D internal noding
!> \param NM Mesh index

SUBROUTINE ADJUST_HT3D_WALL_CELLS(NM)

INTEGER, INTENT(IN) :: NM
INTEGER :: IW,ITW
TYPE(MESH_TYPE), POINTER :: M

M =>  MESHES(NM)

PRIMARY_WALL_LOOP_1: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
   CALL REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL=IW)
ENDDO PRIMARY_WALL_LOOP_1

PRIMARY_THIN_WALL_LOOP_1: DO ITW=1,M%N_THIN_WALL_CELLS
   CALL REALLOCATE_ONE_D_ARRAYS(NM,THIN_WALL_CELL=ITW)
ENDDO PRIMARY_THIN_WALL_LOOP_1

END SUBROUTINE ADJUST_HT3D_WALL_CELLS


!> \brief For a given WALL or THIN_WALL with VARIABLE_THICKNESS or HT3D, adjust the 1-D internal noding
!> \param NM Mesh index
!> \param WALL_CELL Optional WALL cell index
!> \param THIN_WALL_CELL Optional THIN_WALL cell index

SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL)

USE GEOMETRY_FUNCTIONS, ONLY: GET_N_LAYER_CELLS,GET_WALL_NODE_COORDINATES,GET_WALL_NODE_WEIGHTS
USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_REAL_ARRAY,REALLOCATE_INTEGER_ARRAY,PACK_WALL,PACK_THIN_WALL
INTEGER, INTENT(IN) :: NM
INTEGER, INTENT(IN), OPTIONAL :: WALL_CELL,THIN_WALL_CELL
INTEGER :: NL,N_CELLS_MAX,II,NWP,N,I,ITMP,NN
INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX
INTEGER, ALLOCATABLE, DIMENSION(:) :: N_LAYER_CELLS_OLD
REAL(EB) :: MIN_DENSITY
REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_DENSITY
TYPE(MATERIAL_TYPE), POINTER :: ML
REAL(EB), ALLOCATABLE, DIMENSION(:) :: X_S_OLD
LOGICAL, ALLOCATABLE, DIMENSION(:) :: REMESH_LAYER
LOGICAL :: SET_N_CELLS_OLD
TYPE(WALL_TYPE), POINTER :: WC
TYPE(THIN_WALL_TYPE), POINTER :: TW
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(MESH_TYPE), POINTER :: M
TYPE(OBSTRUCTION_TYPE), POINTER :: OB
TYPE(STORAGE_TYPE), POINTER :: OS_DUMMY

M => MESHES(NM)

IF (PRESENT(WALL_CELL)) THEN
   WC => M%WALL(WALL_CELL)
   SF => SURFACE(WC%SURF_INDEX)
   IF (.NOT.SF%VARIABLE_THICKNESS .AND. .NOT.SF%HT_DIM>1) RETURN
   IF (WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) RETURN
   ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX)
   BC  => M%BOUNDARY_COORD(WC%BC_INDEX)
   OB => M%OBSTRUCTION(WC%OBST_INDEX)
ELSEIF (PRESENT(THIN_WALL_CELL)) THEN
   TW => M%THIN_WALL(THIN_WALL_CELL)
   SF => SURFACE(TW%SURF_INDEX)
   IF (SF%HT_DIM==1) RETURN
   ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX)
   BC  => M%BOUNDARY_COORD(TW%BC_INDEX)
   OB => M%OBSTRUCTION(TW%OBST_INDEX)
ENDIF

! This code is probably unnecessary. It is only in case the thickness of the solid has not been determined.

IF (ONE_D%LAYER_THICKNESS(1)<TWO_EPSILON_EB) THEN
   SELECT CASE(ABS(BC%IOR))
      CASE(1) ; ONE_D%LAYER_THICKNESS(1) = OB%X2 - OB%X1
      CASE(2) ; ONE_D%LAYER_THICKNESS(1) = OB%Y2 - OB%Y1
      CASE(3) ; ONE_D%LAYER_THICKNESS(1) = OB%Z2 - OB%Z1
   END SELECT
ENDIF

! Initially these arrays were allocated under the assumption that there is only one layer.

IF (ONE_D%N_LAYERS>1) THEN
   CALL REALLOCATE_REAL_ARRAY(ONE_D%SMALLEST_CELL_SIZE,1,1,ONE_D%N_LAYERS)
   CALL REALLOCATE_REAL_ARRAY(ONE_D%DDSUM,1,1,ONE_D%N_LAYERS)
   CALL REALLOCATE_INTEGER_ARRAY(ONE_D%N_LAYER_CELLS,1,1,ONE_D%N_LAYERS)
   IF (ALLOCATED(ONE_D%REMESH_NWP)) CALL REALLOCATE_INTEGER_ARRAY(ONE_D%REMESH_NWP,1,1,ONE_D%N_LAYERS)
   IF (ALLOCATED(ONE_D%LAYER_THICKNESS_OLD)) CALL REALLOCATE_REAL_ARRAY (ONE_D%LAYER_THICKNESS_OLD,1,1,ONE_D%N_LAYERS)
ENDIF

IF (ALLOCATED(ONE_D%MIN_DIFFUSIVITY)) DEALLOCATE(ONE_D%MIN_DIFFUSIVITY) ; ALLOCATE(ONE_D%MIN_DIFFUSIVITY(1:ONE_D%N_LAYERS))

! Go through all layers and reallocate arrays where necessary

ONE_D%N_CELLS_INI = 0
ONE_D%N_CELLS_MAX = 0
LAYER_DENSITY = 0._EB
SET_N_CELLS_OLD = .FALSE.

LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS

   ! Get the minimum thermal diffusivity for this layer

   ONE_D%MIN_DIFFUSIVITY(NL) = HUGE_EB
   DO NN=1,ONE_D%N_MATL
      IF (ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL)>TWO_EPSILON_EB) &
         ONE_D%MIN_DIFFUSIVITY(NL) = MIN(ONE_D%MIN_DIFFUSIVITY(NL),MATERIAL(ONE_D%MATL_INDEX(NN))%THERMAL_DIFFUSIVITY)
      LAYER_DENSITY(NL) = LAYER_DENSITY(NL) + ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL) / MATERIAL(ONE_D%MATL_INDEX(NN))%RHO_S
   ENDDO
   LAYER_DENSITY(NL) = 1._EB/LAYER_DENSITY(NL)

   ! Get the new N_CELLS_MAX for this wall cell

   MIN_DENSITY = LAYER_DENSITY(NL)
   MATL_REAC_LOOP: DO N=1,ONE_D%N_MATL
      IF (MATERIAL(ONE_D%MATL_INDEX(N))%N_REACTIONS > 0) THEN ! Only compute a swell ratio if there is a reacting material
         SET_N_CELLS_OLD = .TRUE.
         DO NN=1,ONE_D%N_MATL
            IF (MATERIAL(ONE_D%MATL_INDEX(NN))%ALLOW_SWELLING) &
               MIN_DENSITY = MIN(MIN_DENSITY,MATERIAL(ONE_D%MATL_INDEX(NN))%RHO_S)
         ENDDO
         EXIT MATL_REAC_LOOP
      ENDIF
   ENDDO MATL_REAC_LOOP

   CALL GET_N_LAYER_CELLS(ONE_D%MIN_DIFFUSIVITY(NL),LAYER_DENSITY(NL)/MIN_DENSITY*ONE_D%LAYER_THICKNESS(NL), &
                          ONE_D%STRETCH_FACTOR(NL),ONE_D%CELL_SIZE_FACTOR(NL),ONE_D%CELL_SIZE(NL),ONE_D%N_LAYER_CELLS_MAX(NL),&
                          N_CELLS_MAX,ONE_D%SMALLEST_CELL_SIZE(NL),ONE_D%DDSUM(NL))

   ONE_D%N_CELLS_MAX = ONE_D%N_CELLS_MAX + N_CELLS_MAX

   ! Get the new N_CELLS_INI for this wall cell

   CALL GET_N_LAYER_CELLS(ONE_D%MIN_DIFFUSIVITY(NL),ONE_D%LAYER_THICKNESS(NL),ONE_D%STRETCH_FACTOR(NL), &
                          ONE_D%CELL_SIZE_FACTOR(NL),ONE_D%CELL_SIZE(NL),ONE_D%N_LAYER_CELLS_MAX(NL),ONE_D%N_LAYER_CELLS(NL),&
                          ONE_D%SMALLEST_CELL_SIZE(NL),ONE_D%DDSUM(NL))

   ONE_D%N_CELLS_INI = ONE_D%N_CELLS_INI + ONE_D%N_LAYER_CELLS(NL)

ENDDO LAYER_LOOP

NWP_MAX = MAX(NWP_MAX,ONE_D%N_CELLS_MAX)
IF (SET_N_CELLS_OLD) ONE_D%N_CELLS_OLD = ONE_D%N_CELLS_MAX

ALLOCATE(LAYER_INDEX(0:ONE_D%N_CELLS_MAX+1))
NL = 1
DO II=1,ONE_D%N_CELLS_INI
   IF (II>SUM(ONE_D%N_LAYER_CELLS(1:NL))) NL = NL + 1
   LAYER_INDEX(II) = NL
ENDDO
LAYER_INDEX(0) = 1
LAYER_INDEX(ONE_D%N_CELLS_INI+1) = ONE_D%N_LAYERS

IF (ALLOCATED(ONE_D%M_DOT_S_PP)) DEALLOCATE(ONE_D%M_DOT_S_PP) ; ALLOCATE(ONE_D%M_DOT_S_PP(0:ONE_D%N_MATL))
IF (ALLOCATED(ONE_D%X))          DEALLOCATE(ONE_D%X)          ; ALLOCATE(ONE_D%X(0:ONE_D%N_CELLS_MAX))
IF (ALLOCATED(ONE_D%TMP))        DEALLOCATE(ONE_D%TMP)        ; ALLOCATE(ONE_D%TMP(0:ONE_D%N_CELLS_MAX+1))
IF (ALLOCATED(ONE_D%DELTA_TMP))  DEALLOCATE(ONE_D%DELTA_TMP)  ; ALLOCATE(ONE_D%DELTA_TMP(0:ONE_D%N_CELLS_MAX+1))
IF (ALLOCATED(ONE_D%RHO_C_S))    DEALLOCATE(ONE_D%RHO_C_S)    ; ALLOCATE(ONE_D%RHO_C_S(ONE_D%N_CELLS_MAX))
IF (ALLOCATED(ONE_D%K_S))        DEALLOCATE(ONE_D%K_S)        ; ALLOCATE(ONE_D%K_S(0:ONE_D%N_CELLS_MAX+1))
DO NN=1,ONE_D%N_MATL
   IF (ALLOCATED(ONE_D%MATL_COMP(NN)%RHO)) DEALLOCATE(ONE_D%MATL_COMP(NN)%RHO)
   ALLOCATE(ONE_D%MATL_COMP(NN)%RHO(0:ONE_D%N_CELLS_MAX+1))
ENDDO

! Get the new cell coordinates

ALLOCATE(X_S_OLD(0:1)); X_S_OLD=0._EB
ALLOCATE(N_LAYER_CELLS_OLD(ONE_D%N_LAYERS)) ; N_LAYER_CELLS_OLD=1
ALLOCATE(REMESH_LAYER(ONE_D%N_LAYERS)) ; REMESH_LAYER=.TRUE.
CALL GET_WALL_NODE_COORDINATES(ONE_D%N_CELLS_INI,1,ONE_D%N_LAYERS,ONE_D%N_LAYER_CELLS, &
                               N_LAYER_CELLS_OLD,ONE_D%SMALLEST_CELL_SIZE, &
                               ONE_D%STRETCH_FACTOR,REMESH_LAYER,ONE_D%X,X_S_OLD,ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
DEALLOCATE(X_S_OLD)

IF (ALLOCATED(ONE_D%DX_OLD)) DEALLOCATE(ONE_D%DX_OLD) ; ALLOCATE(ONE_D%DX_OLD(ONE_D%N_CELLS_OLD))
DO II=1,MIN(ONE_D%N_CELLS_OLD,ONE_D%N_CELLS_INI)
   ONE_D%DX_OLD(II) = ONE_D%X(II) - ONE_D%X(II-1)
ENDDO

IF (ALLOCATED(ONE_D%LAYER_THICKNESS_OLD)) ONE_D%LAYER_THICKNESS_OLD(1:ONE_D%N_LAYERS) = ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS)

! Reset initial values for some reallocated arrays

ONE_D%TMP = SF%TMP_INNER
ONE_D%DELTA_TMP = 0._EB
ONE_D%K_S = 0._EB

DO II=0,ONE_D%N_CELLS_INI+1
   NL = LAYER_INDEX(II)
   DO NN=1,ONE_D%N_MATL
      ONE_D%MATL_COMP(NN)%RHO(II) = ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL)*LAYER_DENSITY(NL)
   ENDDO
ENDDO

ONE_D%RHO_C_S = 0._EB
NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
POINT_LOOP3: DO I=1,NWP
   ITMP = MIN(I_MAX_TEMP-1,INT(ONE_D%TMP(I)))
   MATERIAL_LOOP3: DO N=1,ONE_D%N_MATL
      IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP3
      ML  => MATERIAL(ONE_D%MATL_INDEX(N))
      ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)
   ENDDO MATERIAL_LOOP3
ENDDO POINT_LOOP3

DEALLOCATE(LAYER_INDEX)

! Count the numbers of REALs, INTEGERs, and LOGICALs in the new WALL or THIN_WALL derived type variable

IF (PRESENT(WALL_CELL)) THEN
   WC%N_REALS=0 ; WC%N_INTEGERS=0 ; WC%N_LOGICALS=0
   CALL PACK_WALL(NM,OS_DUMMY,WC,WC%SURF_INDEX,WC%N_REALS,WC%N_INTEGERS,WC%N_LOGICALS,UNPACK_IT=.FALSE.,COUNT_ONLY=.TRUE.)
ELSEIF (PRESENT(THIN_WALL_CELL)) THEN
   TW%N_INTEGERS=0 ; TW%N_REALS=0
   CALL PACK_THIN_WALL(NM,OS_DUMMY,TW,TW%SURF_INDEX,TW%N_REALS,TW%N_INTEGERS,TW%N_LOGICALS,UNPACK_IT=.FALSE.,COUNT_ONLY=.TRUE.)
ENDIF

END SUBROUTINE REALLOCATE_ONE_D_ARRAYS


!> \brief Set up weighting arrays to transfer 3D solid phase temperatures from one direction sweep to another.
!> \param NM Mesh index

SUBROUTINE INITIALIZE_HT3D_WALL_CELLS(NM)

USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
INTEGER, INTENT(IN) :: NM
INTEGER :: I,IW,IW2,ITW,ITW2,NWP,NWP2,I2,IWA,DM,IOR,NOM,II,JJ,KK,NN,IC,NL
LOGICAL :: IOR_AVOID(-3:3)
REAL(EB) :: X1,X2,Y1,Y2,Z1,Z2,XX1,XX2,YY1,YY2,ZZ1,ZZ2,PRIMARY_VOLUME,OVERLAP_VOLUME,DXX,DYY,DZZ,WEIGHT_FACTOR,&
            SUM_WGT(3),XX,YY,ZZ,WEIGHT,TARGET_WEIGHT
TYPE(WALL_TYPE), POINTER :: WC
TYPE(THIN_WALL_TYPE), POINTER :: TW
TYPE(SURFACE_TYPE), POINTER :: SF,SF2
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D,ONE_D2
TYPE(BOUNDARY_THR_D_TYPE), POINTER :: THR_D
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC,BC2
TYPE(MESH_TYPE), POINTER :: M
INTEGER, ALLOCATABLE, DIMENSION(:) :: INTEGER_DUMMY
REAL(EB), ALLOCATABLE, DIMENSION(:) :: REAL_DUMMY
REAL(EB), PARAMETER :: TOL=0.0001_EB
INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX

M => MESHES(NM)

! Loop over all 3-D wall cells, and for each interior node, find wall or thin wall cells in the other two "alternate"
! coordinate directions

PRIMARY_WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS

   WC => M%WALL(IW)
   SF => SURFACE(WC%SURF_INDEX)

   IF (SF%HT_DIM==1 .OR. WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE PRIMARY_WALL_LOOP

   BC => M%BOUNDARY_COORD(WC%BC_INDEX)
   ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX)
   NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))

   ALLOCATE(LAYER_INDEX(0:ONE_D%N_CELLS_MAX+1))
   NL = 1
   DO II=1,NWP
      IF (II>SUM(ONE_D%N_LAYER_CELLS(1:NL))) NL = NL + 1
      LAYER_INDEX(II) = NL
   ENDDO

   ! Allocate variables that hold information about the wall cells in the two alternate directions

   ALLOCATE(M%BOUNDARY_THR_D(WC%TD_INDEX)%NODE(NWP))
   THR_D => M%BOUNDARY_THR_D(WC%TD_INDEX)
   DO I=1,NWP
      IF (.NOT.ONE_D%HT3D_LAYER(LAYER_INDEX(I))) THEN
         THR_D%NODE(I)%HT3D = .FALSE.
      ELSE
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_INDEX(8))  ; THR_D%NODE(I)%ALTERNATE_WALL_INDEX = 0
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_NODE(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_NODE = 0
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_MESH(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_MESH = 0
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_TYPE(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_TYPE = 0
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_IOR(8))    ; THR_D%NODE(I)%ALTERNATE_WALL_IOR = 0
         ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(8)) ; THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT = 0._EB
      ENDIF
   ENDDO

   DEALLOCATE(LAYER_INDEX)

   X1=BC%X1 ; X2=BC%X2 ; Y1=BC%Y1 ; Y2=BC%Y2 ; Z1=BC%Z1 ; Z2=BC%Z2

   ! Loop over nodes of primary wall cell. For each internal "node", search for
   ! the two alternate wall cells whose 1-D paths intersect the node.

   PRIMARY_NODE_LOOP: DO I=1,NWP

      IF (.NOT.THR_D%NODE(I)%HT3D) CYCLE PRIMARY_NODE_LOOP

      SELECT CASE(BC%IOR)
         CASE( 1) ; X1=BC%X1-ONE_D%X(I)   ; X2=BC%X1-ONE_D%X(I-1)
         CASE(-1) ; X1=BC%X1+ONE_D%X(I-1) ; X2=BC%X1+ONE_D%X(I)
         CASE( 2) ; Y1=BC%Y1-ONE_D%X(I)   ; Y2=BC%Y1-ONE_D%X(I-1)
         CASE(-2) ; Y1=BC%Y1+ONE_D%X(I-1) ; Y2=BC%Y1+ONE_D%X(I)
         CASE( 3) ; Z1=BC%Z1-ONE_D%X(I)   ; Z2=BC%Z1-ONE_D%X(I-1)
         CASE(-3) ; Z1=BC%Z1+ONE_D%X(I-1) ; Z2=BC%Z1+ONE_D%X(I)
      END SELECT
      PRIMARY_VOLUME = (X2-X1)*(Y2-Y1)*(Z2-Z1)
      IOR_AVOID = .FALSE.
      IOR_AVOID(BC%IOR) = .TRUE. ; IOR_AVOID(-BC%IOR) = .TRUE.
      THR_D%NODE(I)%ALTERNATE_WALL_COUNT = 0

      ! Save the mesh number and indices of the mesh cell (II,JJ,KK) in which the center of the solid node is located

      XX = 0.5_EB*(X1+X2) ; YY = 0.5_EB*(Y1+Y2) ; ZZ = 0.5_EB*(Z1+Z2)
      CALL SEARCH_OTHER_MESHES(XX,YY,ZZ,NN,II,JJ,KK)
      IF (NN>0) THEN
         IC = MESHES(NN)%CELL_INDEX(II,JJ,KK)
         IF (MESHES(NN)%CELL(IC)%SOLID) THEN
            THR_D%NODE(I)%I = II
            THR_D%NODE(I)%J = JJ
            THR_D%NODE(I)%K = KK
            THR_D%NODE(I)%MESH_NUMBER = NN
         ENDIF
      ENDIF

      ! Loop over wall cells searching for the "alternate" wall cells whose 1-D path intersects

      ALTERNATE_WALL_LOOP: DO IW2=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS  ! Loop over potential alternate wall cells
         CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NM,WALL_INDEX=IW2)
      ENDDO ALTERNATE_WALL_LOOP

      ALTERNATE_THIN_WALL_LOOP: DO ITW2=1,M%N_THIN_WALL_CELLS  ! Loop over potential alternate wall cells
         CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NM,THIN_WALL_INDEX=ITW2)
      ENDDO ALTERNATE_THIN_WALL_LOOP

      OTHER_MESH_LOOP: DO NOM=1,NMESHES
         IF (NM==NOM) CYCLE
         ALTERNATE_WALL_LOOP_2: DO NN=1,M%OMESH(NOM)%WALL_RECV_BUFFER%N_ITEMS
            IW2 = M%OMESH(NOM)%WALL_RECV_BUFFER%ITEM_INDEX(NN)
            CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,WALL_INDEX=IW2)
         ENDDO ALTERNATE_WALL_LOOP_2
         ALTERNATE_WALL_LOOP_2D: DO NN=1,M%OMESH(NOM)%THIN_WALL_RECV_BUFFER%N_ITEMS  !  THIN_WALL cells, neighboring meshes
            ITW2 = M%OMESH(NOM)%THIN_WALL_RECV_BUFFER%ITEM_INDEX(NN)
            CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,THIN_WALL_INDEX=ITW2)
         ENDDO ALTERNATE_WALL_LOOP_2D
      ENDDO OTHER_MESH_LOOP

      ! Check to see if the HT3D solid object spans the entire width of the computational domain.
      ! There must be at least one exposed surface cell in each coordinate direction.

      DO IOR=1,3
         IF (ABS(BC%IOR)==IOR) CYCLE
         IF (TWO_D .AND. IOR==2) CYCLE
         IF (.NOT.IOR_AVOID(-IOR) .AND. .NOT.IOR_AVOID(IOR)) THEN
            WRITE(LU_ERR,'(7(A,I0))') 'ERROR(423): HT3D solid must have at least one face exposed in direction ',IOR,&
                                      ': Mesh=',NM,', IOR=',BC%IOR,', IIG=',BC%IIG,', JJG=',BC%JJG,', KKG=',BC%KKG,', I=',I
            STOP_STATUS = SETUP_STOP
            RETURN
         ENDIF
      ENDDO

      ! Renormalize weighting factors of the alternate, intersecting 1-D heat conduction paths

      IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT>0 .AND. &
         ABS(SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))-2._EB)>0.001_EB) THEN
         SUM_WGT = 0._EB
         DO IWA=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
            IOR = THR_D%NODE(I)%ALTERNATE_WALL_IOR(IWA)
            SUM_WGT(ABS(IOR)) = SUM_WGT(ABS(IOR)) + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA)
         ENDDO
         DO IWA=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
            IOR = THR_D%NODE(I)%ALTERNATE_WALL_IOR(IWA)
            IF (SUM_WGT(ABS(IOR))>0._EB) THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA) = &
                                         THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA)/SUM_WGT(ABS(IOR))
         ENDDO
         WEIGHT = SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))
         IF (TWO_D) THEN
            TARGET_WEIGHT = 1._EB
         ELSE
            TARGET_WEIGHT = 2._EB
         ENDIF
         IF (ABS(WEIGHT-TARGET_WEIGHT)>0.001_EB) THEN  ! Something is wrong
            WRITE(0,'(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,F6.3)') 'WARNING: Mesh=',NM,' WALL=',IW,' IJK=',BC%IIG,',',&
               BC%JJG,',',BC%KKG,' IOR=',BC%IOR,' NODE=',I,' WEIGHT=',SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))
         ENDIF
      ENDIF

   ENDDO PRIMARY_NODE_LOOP

ENDDO PRIMARY_WALL_LOOP

! Loop over all thin wall cells, all of which are 3-D. When found, find the wall or thin walls in the other coordinate directions.

PRIMARY_THIN_WALL_LOOP: DO ITW=1,M%N_THIN_WALL_CELLS

   TW => M%THIN_WALL(ITW)
   SF => SURFACE(TW%SURF_INDEX)

   IF (SF%HT_DIM==1) CYCLE PRIMARY_THIN_WALL_LOOP

   BC => M%BOUNDARY_COORD(TW%BC_INDEX)
   ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX)
   NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS))
   ALLOCATE(M%BOUNDARY_THR_D(TW%TD_INDEX)%NODE(NWP))
   THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX)
   DO I=1,NWP
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_INDEX(8))  ; THR_D%NODE(I)%ALTERNATE_WALL_INDEX = 0
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_NODE(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_NODE = 0
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_MESH(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_MESH = 0
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_TYPE(8))   ; THR_D%NODE(I)%ALTERNATE_WALL_TYPE = 0
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_IOR(8))    ; THR_D%NODE(I)%ALTERNATE_WALL_IOR = 0
      ALLOCATE(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(8)) ; THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT = 0._EB
   ENDDO

   X1=BC%X1 ; X2=BC%X2 ; Y1=BC%Y1 ; Y2=BC%Y2 ; Z1=BC%Z1 ; Z2=BC%Z2

   ! Loop over the internal nodes of the primary thin wall cell

   PRIMARY_THIN_NODE_LOOP: DO I=1,NWP

      SELECT CASE(BC%IOR)
         CASE( 1) ; X1=BC%X1-ONE_D%X(I)   ; X2=BC%X1-ONE_D%X(I-1)
         CASE(-1) ; X1=BC%X1+ONE_D%X(I-1) ; X2=BC%X1+ONE_D%X(I)
         CASE( 2) ; Y1=BC%Y1-ONE_D%X(I)   ; Y2=BC%Y1-ONE_D%X(I-1)
         CASE(-2) ; Y1=BC%Y1+ONE_D%X(I-1) ; Y2=BC%Y1+ONE_D%X(I)
         CASE( 3) ; Z1=BC%Z1-ONE_D%X(I)   ; Z2=BC%Z1-ONE_D%X(I-1)
         CASE(-3) ; Z1=BC%Z1+ONE_D%X(I-1) ; Z2=BC%Z1+ONE_D%X(I)
      END SELECT
      PRIMARY_VOLUME = (X2-X1)*(Y2-Y1)*(Z2-Z1)
      IOR_AVOID = .FALSE.
      IOR_AVOID(BC%IOR) = .TRUE. ; IOR_AVOID(-BC%IOR) = .TRUE.
      THR_D%NODE(I)%ALTERNATE_WALL_COUNT = 0

      ALTERNATE_WALL_LOOP_B: DO IW2=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS  ! WALL cells, current mesh
         CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NM,WALL_INDEX=IW2)
      ENDDO ALTERNATE_WALL_LOOP_B
      ALTERNATE_THIN_WALL_LOOP_B: DO ITW2=1,M%N_THIN_WALL_CELLS  ! THIN_WALL cells, current mesh
         CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NM,THIN_WALL_INDEX=ITW2)
      ENDDO ALTERNATE_THIN_WALL_LOOP_B

      OTHER_MESH_LOOP_B: DO NOM=1,NMESHES
         IF (NM==NOM) CYCLE
         ALTERNATE_WALL_LOOP_2B: DO NN=1,M%OMESH(NOM)%WALL_RECV_BUFFER%N_ITEMS  ! WALL cells, neighboring meshes
            IW2 = M%OMESH(NOM)%WALL_RECV_BUFFER%ITEM_INDEX(NN)
            CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,WALL_INDEX=IW2)
         ENDDO ALTERNATE_WALL_LOOP_2B
         ALTERNATE_WALL_LOOP_2C: DO NN=1,M%OMESH(NOM)%THIN_WALL_RECV_BUFFER%N_ITEMS  !  THIN_WALL cells, neighboring meshes
            ITW2 = M%OMESH(NOM)%THIN_WALL_RECV_BUFFER%ITEM_INDEX(NN)
            CALL SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,THIN_WALL_INDEX=ITW2)
         ENDDO ALTERNATE_WALL_LOOP_2C
      ENDDO OTHER_MESH_LOOP_B

      ! Check to see if the HT3D solid object spans the entire width of the computational domain.
      ! There must be at least one exposed surface cell in each coordinate direction.

      DO IOR=1,3
         IF (ABS(BC%IOR)==IOR) CYCLE
         IF (.NOT.IOR_AVOID(-IOR) .AND. .NOT.IOR_AVOID(IOR)) THEN
            WRITE(LU_ERR,'(7(A,I0))') 'ERROR(424): HT3D thin solid must have at least one face exposed in direction ',IOR,&
                                      ': Mesh=',NM,', IOR=',BC%IOR,', IIG=',BC%IIG,', JJG=',BC%JJG,', KKG=',BC%KKG,', I=',I
            STOP_STATUS = SETUP_STOP
            RETURN
         ENDIF
      ENDDO

      ! Renormalize the weighting factors for the temperature interpolation

      IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT>0 .AND. &
         ABS(SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))-2._EB)>0.001_EB) THEN
         SUM_WGT = 0._EB
         DO IWA=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
            IOR = THR_D%NODE(I)%ALTERNATE_WALL_IOR(IWA)
            SUM_WGT(ABS(IOR)) = SUM_WGT(ABS(IOR)) + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA)
         ENDDO
         DO IWA=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
            IOR = THR_D%NODE(I)%ALTERNATE_WALL_IOR(IWA)
            IF (SUM_WGT(ABS(IOR))>0._EB) THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA) = &
                                         THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA)/SUM_WGT(ABS(IOR))
         ENDDO
         WEIGHT = SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))
         IF (TWO_D) THEN
            TARGET_WEIGHT = 1._EB
         ELSE
            TARGET_WEIGHT = 2._EB
         ENDIF
         IF (ABS(WEIGHT-TARGET_WEIGHT)>0.001_EB) THEN  ! Something is wrong
            WRITE(0,'(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,F6.3)') 'WARNING: Mesh=',NM,' THIN_WALL=',ITW,' IJK=',BC%IIG,',',&
               BC%JJG,',',BC%KKG,' IOR=',BC%IOR,' NODE=',I,' WEIGHT=',SUM(THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(:))
         ENDIF
      ENDIF

   ENDDO PRIMARY_THIN_NODE_LOOP

ENDDO PRIMARY_THIN_WALL_LOOP

CONTAINS

!> \brief Find WALL or THIN_WALL cells whose internal nodes overlap those of the primary WALL or THIN_WALL cell
!> \param NOM Mesh number of the primary cell or its neighbor
!> \param WALL_INDEX Optional wall cell index
!> \param THIN_WALL_INDEX Optional thin wall cell index

SUBROUTINE SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,WALL_INDEX,THIN_WALL_INDEX)

TYPE (MESH_TYPE), POINTER :: M2
INTEGER, INTENT(IN) :: NOM
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,THIN_WALL_INDEX
INTEGER :: CELL
TYPE(WALL_TYPE), POINTER :: WC2
TYPE(THIN_WALL_TYPE), POINTER :: TW2

M2  => MESHES(NOM)

IF (PRESENT(WALL_INDEX)) THEN
   CELL = WALL_INDEX
   WC2 => M2%WALL(WALL_INDEX)
   SF2 => SURFACE(WC2%SURF_INDEX)
   IF (SF2%HT_DIM==1 .OR. WC2%BOUNDARY_TYPE/=SOLID_BOUNDARY) RETURN
   BC2 => M2%BOUNDARY_COORD(WC2%BC_INDEX)
   ONE_D2 => M2%BOUNDARY_ONE_D(WC2%OD_INDEX)
ELSE
   CELL = THIN_WALL_INDEX
   TW2 => M2%THIN_WALL(THIN_WALL_INDEX)
   SF2 => SURFACE(TW2%SURF_INDEX)
   BC2 => M2%BOUNDARY_COORD(TW2%BC_INDEX)
   ONE_D2 => M2%BOUNDARY_ONE_D(TW2%OD_INDEX)
ENDIF

IF (IOR_AVOID(BC2%IOR)) RETURN
IF (PRESENT(WALL_INDEX) .AND. (BC2%IIG==0 .OR. BC2%IIG==M2%IBP1)) RETURN
IF (PRESENT(WALL_INDEX) .AND. (BC2%JJG==0 .OR. BC2%JJG==M2%JBP1)) RETURN
IF (PRESENT(WALL_INDEX) .AND. (BC2%KKG==0 .OR. BC2%KKG==M2%KBP1)) RETURN

XX1=BC2%X1 ; XX2=BC2%X2 ; YY1=BC2%Y1 ; YY2=BC2%Y2 ; ZZ1=BC2%Z1 ; ZZ2=BC2%Z2

IF (ABS(BC2%IOR)/=1) THEN ; DXX = MIN(XX2,X2)-MAX(XX1,X1) ; IF (DXX<=0._EB) RETURN ; ENDIF
IF (ABS(BC2%IOR)/=2) THEN ; DYY = MIN(YY2,Y2)-MAX(YY1,Y1) ; IF (DYY<=0._EB) RETURN ; ENDIF
IF (ABS(BC2%IOR)/=3) THEN ; DZZ = MIN(ZZ2,Z2)-MAX(ZZ1,Z1) ; IF (DZZ<=0._EB) RETURN ; ENDIF

NWP2 = SUM(ONE_D2%N_LAYER_CELLS(1:ONE_D2%N_LAYERS))

ALTERNATE_NODE_LOOP: DO I2=1,NWP2  ! Loop over nodes of alternate wall cell
   SELECT CASE(BC2%IOR)
      CASE( 1) ; XX1=BC2%X2-ONE_D2%X(I2)   ; XX2=BC2%X2-ONE_D2%X(I2-1)
                 DXX = MIN(XX2,X2)-MAX(XX1,X1) ; IF (DXX<=0._EB) CYCLE ALTERNATE_NODE_LOOP
      CASE(-1) ; XX1=BC2%X1+ONE_D2%X(I2-1) ; XX2=BC2%X1+ONE_D2%X(I2)
                 DXX = MIN(XX2,X2)-MAX(XX1,X1) ; IF (DXX<=0._EB) CYCLE ALTERNATE_NODE_LOOP
      CASE( 2) ; YY1=BC2%Y2-ONE_D2%X(I2)   ; YY2=BC2%Y2-ONE_D2%X(I2-1)
                 DYY = MIN(YY2,Y2)-MAX(YY1,Y1) ; IF (DYY<=0._EB) CYCLE ALTERNATE_NODE_LOOP
      CASE(-2) ; YY1=BC2%Y1+ONE_D2%X(I2-1) ; YY2=BC2%Y1+ONE_D2%X(I2)
                 DYY = MIN(YY2,Y2)-MAX(YY1,Y1) ; IF (DYY<=0._EB) CYCLE ALTERNATE_NODE_LOOP
      CASE( 3) ; ZZ1=BC2%Z2-ONE_D2%X(I2)   ; ZZ2=BC2%Z2-ONE_D2%X(I2-1)
                 DZZ = MIN(ZZ2,Z2)-MAX(ZZ1,Z1) ; IF (DZZ<=0._EB) CYCLE ALTERNATE_NODE_LOOP
      CASE(-3) ; ZZ1=BC2%Z1+ONE_D2%X(I2-1) ; ZZ2=BC2%Z1+ONE_D2%X(I2)
                 DZZ = MIN(ZZ2,Z2)-MAX(ZZ1,Z1) ; IF (DZZ<=0._EB) CYCLE ALTERNATE_NODE_LOOP
   END SELECT
   OVERLAP_VOLUME = DXX*DYY*DZZ
   WEIGHT_FACTOR = OVERLAP_VOLUME/PRIMARY_VOLUME
   IF (WEIGHT_FACTOR<TOL) CYCLE ALTERNATE_NODE_LOOP
   DM = SIZE(THR_D%NODE(I)%ALTERNATE_WALL_INDEX)
   IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT+1>DM) CALL REALLOCATE_ALTERNATE
   THR_D%NODE(I)%ALTERNATE_WALL_COUNT = THR_D%NODE(I)%ALTERNATE_WALL_COUNT + 1
   IWA = THR_D%NODE(I)%ALTERNATE_WALL_COUNT
   THR_D%NODE(I)%ALTERNATE_WALL_MESH(IWA) = NOM
   THR_D%NODE(I)%ALTERNATE_WALL_INDEX(IWA) = CELL
   THR_D%NODE(I)%ALTERNATE_WALL_IOR(IWA) = BC2%IOR
   IF (PRESENT(THIN_WALL_INDEX)) THR_D%NODE(I)%ALTERNATE_WALL_TYPE(IWA) = 1
   THR_D%NODE(I)%ALTERNATE_WALL_NODE(IWA) = I2
   THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(IWA) = WEIGHT_FACTOR
   IOR_AVOID(-BC2%IOR) = .TRUE.  ! Do not use the opposite side wall
ENDDO ALTERNATE_NODE_LOOP

END SUBROUTINE SEARCH_FOR_ALTERNATE_WALL_CELLS


SUBROUTINE REALLOCATE_ALTERNATE

ALLOCATE(INTEGER_DUMMY(DM+8))
INTEGER_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_MESH(1:DM)
CALL MOVE_ALLOC(INTEGER_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_MESH)
THR_D%NODE(I)%ALTERNATE_WALL_MESH(DM+1:DM+8) = 0
ALLOCATE(INTEGER_DUMMY(DM+8))
INTEGER_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_IOR(1:DM)
CALL MOVE_ALLOC(INTEGER_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_IOR)
THR_D%NODE(I)%ALTERNATE_WALL_IOR(DM+1:DM+8) = 0
ALLOCATE(INTEGER_DUMMY(DM+8))
INTEGER_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_TYPE(1:DM)
CALL MOVE_ALLOC(INTEGER_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_TYPE)
THR_D%NODE(I)%ALTERNATE_WALL_TYPE(DM+1:DM+8) = 0
ALLOCATE(INTEGER_DUMMY(DM+8))
INTEGER_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(1:DM)
CALL MOVE_ALLOC(INTEGER_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_INDEX)
THR_D%NODE(I)%ALTERNATE_WALL_INDEX(DM+1:DM+8) = 0
ALLOCATE(INTEGER_DUMMY(DM+8))
INTEGER_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_NODE(1:DM)
CALL MOVE_ALLOC(INTEGER_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_NODE)
ALLOCATE(REAL_DUMMY(DM+8))
REAL_DUMMY(1:DM) = THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(1:DM)
CALL MOVE_ALLOC(REAL_DUMMY,THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT)
THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(DM+1:DM+8)=0._EB

END SUBROUTINE REALLOCATE_ALTERNATE

END SUBROUTINE INITIALIZE_HT3D_WALL_CELLS


!> \brief Initialize a few GEOM arrays

SUBROUTINE INITIALIZE_MESH_VARIABLES_3(NM)

INTEGER, INTENT(IN) :: NM
INTEGER :: ICF,N_REALS,N_INTEGERS,N_LOGICALS
TYPE(CFACE_TYPE), POINTER :: CFA
TYPE(MESH_TYPE), POINTER :: M

M => MESHES(NM)

IF (M%N_CFACE_CELLS_DIM<1) RETURN

N_REALS    = 0
N_INTEGERS = 0
N_LOGICALS = 0
DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTWALL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
   CFA => M%CFACE(ICF)
   N_REALS    = MAX(N_REALS   ,CFA%N_REALS)
   N_INTEGERS = MAX(N_INTEGERS,CFA%N_INTEGERS)
   N_LOGICALS = MAX(N_LOGICALS,CFA%N_LOGICALS)
ENDDO

ALLOCATE(M%CFACE_STORAGE%REALS(N_REALS))
ALLOCATE(M%CFACE_STORAGE%INTEGERS(N_INTEGERS))
ALLOCATE(M%CFACE_STORAGE%LOGICALS(N_LOGICALS))

END SUBROUTINE INITIALIZE_MESH_VARIABLES_3


!> \brief Intialize Crayfishpak (FFT) Poisson solver
!> \param NM Mesh number

SUBROUTINE INITIALIZE_POISSON_SOLVER(NM)

USE POIS, ONLY: H3CZIS,H2CZIS,H3CSIS,H2CYIS
INTEGER, INTENT(IN) :: NM
REAL(EB) :: XLM,XMU,XS,YS,ZS,XF,YF,ZF
INTEGER  :: N,IZERO,IERR,IBP1,JBP1,KBP1,IBAR,JBAR,KBAR,IW,IOR,JDIM
INTEGER, POINTER :: ITRN,JTRN,KTRN,LBC,MBC,NBC
INTEGER, POINTER, DIMENSION(:) :: NOC
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (WALL_TYPE), POINTER :: WC
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (MESH_TYPE), POINTER :: M
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC

IERR = 0
M => MESHES(NM)
IBP1 = M%IBP1
JBP1 = M%JBP1
KBP1 = M%KBP1
IBAR = M%IBAR
JBAR = M%JBAR
KBAR = M%KBAR
XS = M%XS
YS = M%YS
ZS = M%ZS
XF = M%XF
YF = M%YF
ZF = M%ZF

! Allocate major arrays

ITRN =>M%ITRN
JTRN =>M%JTRN
KTRN =>M%KTRN
LBC =>M%LBC
MBC =>M%MBC
NBC =>M%NBC
NOC=>TRANS(NM)%NOC
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=0
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=1
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=2
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=3
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=4
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=5
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) M%IPS=6
SELECT CASE(PRES_FLAG)
   CASE DEFAULT
      IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) THEN
         WRITE(LU_ERR,'(A,I0,A)') 'ERROR(425): MESH ',NM,' can stretch in at most 2 coordinate directions.'
         STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF
   CASE (GLMAT_FLAG,UGLMAT_FLAG)
      M%IPS=0 ! For ULMAT_FLAG, IPS set to 0 (no transpose of indices) for ZM%USE_FFT=F in ULMAT_SOLVER_SETUP
END SELECT

IF (M%IPS<=1 .OR. M%IPS==4) THEN
   ITRN = IBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = KBP1

   ! pressure periodic boundary conditions
   IF (FISHPAK_BC(1)==FISHPAK_BC_PERIODIC) ITRN=IBAR
   IF (FISHPAK_BC(2)==FISHPAK_BC_PERIODIC) JTRN=JBAR
   IF (FISHPAK_BC(3)==FISHPAK_BC_PERIODIC) KTRN=KBAR
ENDIF

IF (M%IPS==2) THEN
   ITRN = JBP1
   JTRN = IBP1
   KTRN = KBP1
   ALLOCATE(M%BZST(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF

IF (M%IPS==3 .OR. M%IPS==6) THEN
   ITRN = KBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = IBP1
   ALLOCATE(M%BXST(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
   ALLOCATE(M%BYST(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYST',IZERO)
   ALLOCATE(M%BYFT(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYFT',IZERO)
   ALLOCATE(M%BZST(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF

IF (M%IPS==5) THEN
   ITRN = IBP1
   JTRN = KBP1
   KTRN = JBP1
   ALLOCATE(M%BXST(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
ENDIF

IF (M%IPS==7) THEN
   ITRN = IBP1
   JTRN = JBP1
   KTRN = 1
ENDIF

IF (M%IPS<=3 .OR. M%IPS==7) THEN
   M%LSAVE = (ITRN+1)*JTRN*KTRN+7*ITRN+5*JTRN+6*KTRN+56
   M%LWORK = (ITRN+1)*JTRN*KTRN
ELSE
   N_LOOP: DO N=1,50
      IF ((JTRN+1)<=2**N) EXIT N_LOOP
   ENDDO N_LOOP
   M%LSAVE = KTRN*(6*N*(2**N)+2*N+19)+8*ITRN+7*JTRN+38
   M%LWORK = JTRN*(ITRN*(KTRN+1)+1)
ENDIF

ALLOCATE(M%SAVE1(-3:M%LSAVE),STAT=IZERO)    ; CALL ChkMemErr('INIT','SAVE1',IZERO)
ALLOCATE(M%WORK(M%LWORK),STAT=IZERO)        ; CALL ChkMemErr('INIT','WORK',IZERO)
ALLOCATE(M%PRHS(ITRN,JTRN,KTRN),STAT=IZERO) ; CALL ChkMemErr('INIT','PRHS',IZERO)
IF (JBAR>1 ) JDIM = JBP1
IF (JBAR==1) JDIM = 1
ALLOCATE(M%BXS(JDIM,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','BXS',IZERO)
ALLOCATE(M%BXF(JDIM,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','BXF',IZERO)
ALLOCATE(M%BYS(IBP1,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','BYS',IZERO)
ALLOCATE(M%BYF(IBP1,KBP1),STAT=IZERO) ; CALL ChkMemErr('INIT','BYF',IZERO)
ALLOCATE(M%BZS(IBP1,JDIM),STAT=IZERO) ; CALL ChkMemErr('INIT','BZS',IZERO)
ALLOCATE(M%BZF(IBP1,JDIM),STAT=IZERO) ; CALL ChkMemErr('INIT','BZF',IZERO)

M%POIS_PTB = 0._EB
M%SAVE1 = 0._EB
M%WORK = 0._EB
M%PRHS = 0._EB
M%BXS  = 0._EB
M%BXF  = 0._EB
M%BYS  = 0._EB
M%BYF  = 0._EB
M%BZS  = 0._EB
M%BZF  = 0._EB

! Initialize pressure solver

XLM = 0._EB         ! No Helmholtz equation
XMU = 0._EB         ! No Helmholtz equation

! BC index for Fishpak solver

! From Roland Sweet's notes:
!
! Here we use LBC as an example, this is the BC index for the X direction. MBC and NMC are
! analogous for the Y and Z directions.
!
! LBC = 0, solution is periodic in X.
!
! LBC = 1, solution is specified at XS (Dirichlet) and XF (Dirichlet).
!
! LBC = 2, solution is specified at XS (Dirichlet) and derivative of solution is specified at XF (Neumann).
!
! LBC = 3, derivative of solution is specified at XS (Neumann) and XF (Neumann).
!
! LBC = 4, derivative of solution is specified at XS (Neumann) and solution is specified at XF (Dirichlet).
!
! LBC = 5, the solution is unspecified at r = RS = 0 and the solution is specified at r = RF.
!
! LBC = 6, if the solution is unspecified at r = RS = 0 and the derivative of the solution with respect to r is specified
!          at r = RF.

LBC = FISHPAK_BC_NEUMANN_NEUMANN
MBC = FISHPAK_BC_NEUMANN_NEUMANN
NBC = FISHPAK_BC_NEUMANN_NEUMANN

! Look for OPEN vents -- this will change the entire face to DIRICHLET BCs

VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   IF (VT%BOUNDARY_TYPE /= OPEN_BOUNDARY) CYCLE VENT_LOOP
   IF (VT%I1==0 .AND. VT%I2==0) THEN
      IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (LBC==FISHPAK_BC_NEUMANN_DIRICHLET) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%I1==M%IBAR .AND. VT%I2==M%IBAR) THEN
      IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%J1==0 .AND. VT%J2==0) THEN
      IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (MBC==FISHPAK_BC_NEUMANN_DIRICHLET) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%J1==M%JBAR .AND. VT%J2==M%JBAR) THEN
      IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (MBC==FISHPAK_BC_DIRICHLET_NEUMANN) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%K1==0 .AND. VT%K2==0) THEN
      IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_DIRICHLET_NEUMANN
      IF (NBC==FISHPAK_BC_NEUMANN_DIRICHLET) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
   IF (VT%K1==M%KBAR .AND. VT%K2==M%KBAR) THEN
      IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_NEUMANN_DIRICHLET
      IF (NBC==FISHPAK_BC_DIRICHLET_NEUMANN) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   ENDIF
ENDDO VENT_LOOP

! All interpolated boundaries are Dirichlet

DO IW=1,M%N_EXTERNAL_WALL_CELLS
   WC => M%WALL(IW)
   BC => M%BOUNDARY_COORD(WC%BC_INDEX)
   IF (M%EXTERNAL_WALL(IW)%NOM==0) CYCLE
   SELECT CASE(BC%IOR)
      CASE( 1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (LBC==FISHPAK_BC_NEUMANN_DIRICHLET) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN)   LBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN) LBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE( 2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (MBC==FISHPAK_BC_NEUMANN_DIRICHLET) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN)   MBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (MBC==FISHPAK_BC_DIRICHLET_NEUMANN) MBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE( 3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_DIRICHLET_NEUMANN
         IF (NBC==FISHPAK_BC_NEUMANN_DIRICHLET) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
      CASE(-3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN)   NBC = FISHPAK_BC_NEUMANN_DIRICHLET
         IF (NBC==FISHPAK_BC_DIRICHLET_NEUMANN) NBC = FISHPAK_BC_DIRICHLET_DIRICHLET
   END SELECT
ENDDO

! Periodic pressure boundary conditions for CrayFishpak

IF (FISHPAK_BC(1)==FISHPAK_BC_PERIODIC) LBC=FISHPAK_BC_PERIODIC
IF (FISHPAK_BC(2)==FISHPAK_BC_PERIODIC) MBC=FISHPAK_BC_PERIODIC
IF (FISHPAK_BC(3)==FISHPAK_BC_PERIODIC) NBC=FISHPAK_BC_PERIODIC

! Poisson solver with stretching in the 1st coordinate

SELECT_POISSON_SOLVER: SELECT CASE(M%IPS)

   CASE (0:1) SELECT_POISSON_SOLVER
      IF (.NOT.TWO_D) CALL H3CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,JTRN,IERR,M%SAVE1)
      IF (TWO_D .AND. .NOT.CYLINDRICAL) CALL H2CZIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HX,XLM,ITRN,IERR,M%SAVE1)
      IF (TWO_D .AND. CYLINDRICAL) THEN
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_DIRICHLET_DIRICHLET) LBC = 5
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_DIRICHLET_NEUMANN)   LBC = 6
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_NEUMANN_NEUMANN)     LBC = 6
         IF (ABS(XS)<=TWO_EPSILON_EB .AND. LBC==FISHPAK_BC_NEUMANN_DIRICHLET)   LBC = 5
         CALL H2CYIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,XLM,XMU,ITRN,IERR,M%SAVE1)
      ENDIF
   CASE (2) SELECT_POISSON_SOLVER
      CALL H3CZIS(YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HY,XLM,ITRN,JTRN,IERR,M%SAVE1)
   CASE (3) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CZIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,JTRN,IERR,M%SAVE1)
      ENDIF
   CASE (4) SELECT_POISSON_SOLVER
      CALL H3CSIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HY)
   CASE (5) SELECT_POISSON_SOLVER
      IF (TWO_D) THEN
         CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE1)
      ELSE
         CALL H3CSIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HX,M%HZ)
      ENDIF
   CASE (6) SELECT_POISSON_SOLVER
      CALL H3CSIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,XLM,ITRN,JTRN,IERR,M%SAVE1,M%WORK,M%HZ,M%HY)
   CASE (7) SELECT_POISSON_SOLVER
      CALL H2CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,M%HX,XLM,ITRN,IERR,M%SAVE1)

END SELECT SELECT_POISSON_SOLVER

! Specify the pressure boundary condition for each wall cell

WALL_CELL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS
   WC => M%WALL(IW)
   EWC => M%EXTERNAL_WALL(IW)
   BC => M%BOUNDARY_COORD(WC%BC_INDEX)
   IOR = BC%IOR
   SELECT CASE(IOR)
      CASE( 1)
         IF (LBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. &
             LBC==FISHPAK_BC_NEUMANN_DIRICHLET   .OR. LBC==6) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (LBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
             LBC==FISHPAK_BC_DIRICHLET_NEUMANN   .OR. LBC==5) EWC%PRESSURE_BC_TYPE = DIRICHLET
      CASE(-1)
         IF (LBC==FISHPAK_BC_DIRICHLET_NEUMANN   .OR. &
             LBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. LBC==6) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (LBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. &
             LBC==FISHPAK_BC_NEUMANN_DIRICHLET   .OR. LBC==5) EWC%PRESSURE_BC_TYPE = DIRICHLET
      CASE( 2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. MBC==FISHPAK_BC_NEUMANN_DIRICHLET) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (MBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC==FISHPAK_BC_DIRICHLET_NEUMANN) EWC%PRESSURE_BC_TYPE = DIRICHLET
      CASE(-2)
         IF (MBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. MBC==FISHPAK_BC_DIRICHLET_NEUMANN) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (MBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. MBC==FISHPAK_BC_NEUMANN_DIRICHLET) EWC%PRESSURE_BC_TYPE = DIRICHLET
      CASE( 3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. NBC==FISHPAK_BC_NEUMANN_DIRICHLET) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (NBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC==FISHPAK_BC_DIRICHLET_NEUMANN) EWC%PRESSURE_BC_TYPE = DIRICHLET
      CASE(-3)
         IF (NBC==FISHPAK_BC_NEUMANN_NEUMANN     .OR. NBC==FISHPAK_BC_DIRICHLET_NEUMANN) EWC%PRESSURE_BC_TYPE = NEUMANN
         IF (NBC==FISHPAK_BC_DIRICHLET_DIRICHLET .OR. NBC==FISHPAK_BC_NEUMANN_DIRICHLET) EWC%PRESSURE_BC_TYPE = DIRICHLET
   END SELECT
ENDDO WALL_CELL_LOOP

! Check for errors with Poisson solver initialization

IF (IERR/=0) THEN
   WRITE(LU_ERR,'(A,I0,A,I0)') 'ERROR(426): MESH ',NM,' Poisson initialization error: ',IERR
   STOP_STATUS = SETUP_STOP
   RETURN
ENDIF

END SUBROUTINE INITIALIZE_POISSON_SOLVER


!> \brief Find the WALL_INDEX for a device that is near a solid wall
!> \param NM Mesh number

SUBROUTINE INITIALIZE_DEVICES(NM)

USE COMPLEX_GEOMETRY, ONLY : GET_CFACE_INDEX
INTEGER, INTENT(IN) :: NM
INTEGER :: III,N,II,JJ,KK,IOR,IW,SURF_INDEX,IIG,JJG,KKG,ICF
REAL(EB) :: DEPTH
TYPE (DEVICE_TYPE), POINTER :: DV
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)

DEVICE_LOOP: DO N=1,N_DEVC

   DV => DEVICE(N)

   IF (DV%QUANTITY_INDEX(1)>=0) CYCLE DEVICE_LOOP  ! Do not process gas phsae devices

   IF (DV%INIT_ID=='null' .AND. DV%LP_TAG==0) THEN ! Assume the device is tied to a WALL cell or CFACE

      IF (NM/=DV%MESH) CYCLE DEVICE_LOOP
      II  = INT(GINV(DV%X-M%XS,1,NM)*M%RDXI   + 1._EB)
      JJ  = INT(GINV(DV%Y-M%YS,2,NM)*M%RDETA  + 1._EB)
      KK  = INT(GINV(DV%Z-M%ZS,3,NM)*M%RDZETA + 1._EB)
      IIG = II
      JJG = JJ
      KKG = KK
      IOR = DV%IOR
      IW  = 0
      ICF = 0

      IF (IOR/=0) CALL GET_WALL_INDEX(NM,IIG,JJG,KKG,IOR,IW)

      IF (IW==0 .AND. CC_IBM)  CALL GET_CFACE_INDEX(NM,IIG,JJG,KKG,DV%X,DV%Y,DV%Z,ICF)

      IF (IW==0 .AND. ICF==0 .AND. DV%SPATIAL_STATISTIC=='null') THEN
         WRITE(LU_ERR,'(A,A,A)') 'ERROR(427): DEVC ',TRIM(DV%ID),' requires repositioning.'
         STOP_STATUS = SETUP_STOP
         RETURN
      ELSEIF (IW>0) THEN
         DV%WALL_INDEX = IW
         SURF_INDEX = M%WALL(IW)%SURF_INDEX
      ELSEIF (ICF>0) THEN
         DV%CFACE_INDEX = ICF
         SURF_INDEX = M%CFACE(ICF)%SURF_INDEX
      ELSE
         SURF_INDEX = DV%SURF_INDEX
      ENDIF

   ELSE ! Assume the device is tied to a particle

      IF (DV%PART_CLASS_INDEX<1) CYCLE DEVICE_LOOP
      SURF_INDEX = LAGRANGIAN_PARTICLE_CLASS(DV%PART_CLASS_INDEX)%SURF_INDEX

   ENDIF

   ! Make sure that thermally-thick output is appropriate

   IF (OUTPUT_QUANTITY(DV%QUANTITY_INDEX(1))%INSIDE_SOLID) THEN
      IF (SURFACE(SURF_INDEX)%THERMAL_BC_INDEX /= THERMALLY_THICK) THEN
         WRITE(LU_ERR,'(A,A,A)') 'ERROR(428): DEVC ',TRIM(DV%ID),' must be associated with a heat-conducting surface.'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
      IF (DV%DEPTH>TWO_EPSILON_EB) THEN
         DEPTH = DV%DEPTH
      ELSE
         DEPTH = MAX(0._EB,SUM(SURFACE(SURF_INDEX)%LAYER_THICKNESS)+DV%DEPTH)
      ENDIF
      DV%I_DEPTH = SURFACE(SURF_INDEX)%N_CELLS_INI
      DO III=SURFACE(SURF_INDEX)%N_CELLS_INI,1,-1
         IF (DEPTH<=SURFACE(SURF_INDEX)%X_S(III)) DV%I_DEPTH = III
      ENDDO
   ENDIF

ENDDO DEVICE_LOOP

END SUBROUTINE INITIALIZE_DEVICES


!> \brief Initialize output PROFiles
!> \param NM Mesh number

SUBROUTINE INITIALIZE_PROFILES(NM)

INTEGER, INTENT(IN) :: NM
INTEGER :: NN,N,II,JJ,KK,IW,IOR
LOGICAL :: SUCCESS
TYPE (PROFILE_TYPE), POINTER :: PF
TYPE (MESH_TYPE), POINTER :: M
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
CHARACTER(LABEL_LENGTH) :: HEADING

PROF_LOOP: DO N=1,N_PROF

   PF => PROFILE(N)

   IF (PF%IOR/=0) THEN  ! The PROFile is for a WALL cell
      IF (PF%MESH/=NM) CYCLE PROF_LOOP
      M => MESHES(NM)
      IOR = PF%IOR
      II  = INT(GINV(PF%X-M%XS,1,NM)*M%RDXI   + 1._EB)
      JJ  = INT(GINV(PF%Y-M%YS,2,NM)*M%RDETA  + 1._EB)
      KK  = INT(GINV(PF%Z-M%ZS,3,NM)*M%RDZETA + 1._EB)
      CALL GET_WALL_INDEX(NM,II,JJ,KK,IOR,IW)
      IF (IW>0) THEN
         PF%WALL_INDEX = IW
         SF => SURFACE(M%WALL(IW)%SURF_INDEX)
         ONE_D => M%BOUNDARY_ONE_D(M%WALL(IW)%OD_INDEX)
         BC => M%BOUNDARY_COORD(M%WALL(IW)%BC_INDEX)
      ELSE
         WRITE(LU_ERR,'(A,I0,A)') 'ERROR(429): PROF ',PF%ORDINAL,' requires repositioning.'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
   ELSE  ! The PROFile is for a Lagrangian PARTicle
      SF => SURFACE(LAGRANGIAN_PARTICLE_CLASS(PF%PART_CLASS_INDEX)%SURF_INDEX)
   ENDIF

   ! Check for potential errors

   IF (SF%THERMAL_BC_INDEX/=THERMALLY_THICK) THEN
      WRITE(LU_ERR,'(A,I0,A)') 'ERROR(430): PROF ',N,' must be associated with a heat-conducting surface.'
      STOP_STATUS = SETUP_STOP
      RETURN
   ENDIF

   IF (PF%MATL_INDEX>0) THEN
      SUCCESS = .FALSE.
      DO NN=1,SF%N_MATL
         IF (PF%MATL_INDEX==SF%MATL_INDEX(NN)) THEN
            SUCCESS = .TRUE.
            EXIT
         ENDIF
      ENDDO
      IF (.NOT. SUCCESS) THEN
         WRITE(LU_ERR,'(A,I3,5A)') 'ERROR PROF ',N,'. MATL_ID ',TRIM(MATERIAL(PF%MATL_INDEX)%ID),&
                              ' not part of surface type ',TRIM(SF%ID),' at the profile location.'
         STOP_STATUS = SETUP_STOP
         RETURN
      ENDIF
   ENDIF

   ! If the PROFile is applied to a particle, let the root MPI process open and close the file. Other MPI processes can then
   ! open and write to the file if the particle moves from mesh to mesh.

   IF (PF%IOR==0 .AND. NM>1) CYCLE PROF_LOOP

   IF (APPEND .AND. PF%FORMAT_INDEX==1) THEN
      OPEN(LU_PROF(N),FILE=FN_PROF(N),FORM='FORMATTED',STATUS='OLD',POSITION='APPEND')
   ELSE
      OPEN(LU_PROF(N),FILE=FN_PROF(N),FORM='FORMATTED',STATUS='REPLACE')
      IF (PF%FORMAT_INDEX==1) THEN
         IF (PF%IOR/=0) THEN ! Wall cell
            WRITE(LU_PROF(N),'(A)') "ID, IOR, face center x(m), face center y(m), face center z(m)"
            WRITE(LU_PROF(N),'(A,A,I3,A,E16.9,A,E16.9,A,E16.9)') TRIM(PF%ID),", ",PF%IOR,", ",BC%X,", ",BC%Y,", ",BC%Z
         ELSE
            WRITE(LU_PROF(N),'(A)') TRIM(PF%ID)
         ENDIF
         IF (PF%ID/='null') THEN
            HEADING = PF%ID
         ELSE
            HEADING = OUTPUT_QUANTITY(PF%QUANTITY_INDEX)%SHORT_NAME
         ENDIF
         WRITE(LU_PROF(N),'(A,A)') "Time(s), Npoints, Npoints x Depth (m), Npoints x ",TRIM(HEADING)
      ENDIF
   ENDIF

   CLOSE(LU_PROF(N))

ENDDO PROF_LOOP

END SUBROUTINE INITIALIZE_PROFILES


!> \brief Find the wall index corresponding to the -IOR face of cell (II,JJ,KK)
!> \param NM Mesh number
!> \param II x-index on the inside of the WALL face
!> \param JJ y-index on the inside of the WALL face
!> \param KK z-index on the inside of the WALL face
!> \param IOR Orientation index of the WALL face
!> \param IW Index of the WALL face

SUBROUTINE GET_WALL_INDEX(NM,II,JJ,KK,IOR,IW)

INTEGER, INTENT(IN) :: NM,IOR
INTEGER, INTENT(OUT) :: IW
INTEGER :: IC,II,JJ,KK
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
IC  = M%CELL_INDEX(II,JJ,KK)

IF (M%CELL(IC)%SOLID) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      II = II-1
      CASE( 1)
         IF (II<M%IBP1) II = II+1
      CASE(-2)
         IF (JJ>0)      JJ = JJ-1
      CASE( 2)
         IF (JJ<M%JBP1) JJ = JJ+1
      CASE(-3)
         IF (KK>0)      KK = KK-1
      CASE( 3)
         IF (KK<M%KBP1) KK = KK+1
   END SELECT
ENDIF

IC  = M%CELL_INDEX(II,JJ,KK)
IW  = M%CELL(IC)%WALL_INDEX(-IOR)

IF (IW<=0) THEN
   SELECT CASE(IOR)
      CASE(-1)
         IF (II>0)      IC = M%CELL_INDEX(II-1,JJ,KK)
      CASE( 1)
         IF (II<M%IBP1) IC = M%CELL_INDEX(II+1,JJ,KK)
      CASE(-2)
         IF (JJ>0)      IC = M%CELL_INDEX(II,JJ-1,KK)
      CASE( 2)
         IF (JJ<M%JBP1) IC = M%CELL_INDEX(II,JJ+1,KK)
      CASE(-3)
         IF (KK>0)      IC = M%CELL_INDEX(II,JJ,KK-1)
      CASE( 3)
         IF (KK<M%KBP1) IC = M%CELL_INDEX(II,JJ,KK+1)
   END SELECT
   IW = M%CELL(IC)%WALL_INDEX(-IOR)
ENDIF

END SUBROUTINE GET_WALL_INDEX


!> \brief Initialize time, printout and plot clocks

SUBROUTINE INITIALIZE_GLOBAL_VARIABLES

INTEGER :: IZERO, IG

ICYC              = 0
T_LAST_DUMP_HRR   = T_BEGIN
T_LAST_DUMP_MASS  = T_BEGIN
T_LAST_DUMP_MOM   = T_BEGIN

! N_FACE manages the geometry output time GEOM_CLOCK:

DO IG=1,N_GEOMETRY; N_FACE = N_FACE + GEOMETRY(IG)%N_FACES; ENDDO

ALLOCATE(ENTHALPY_SUM(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','ENTHALPY_SUM',IZERO)
ENTHALPY_SUM = 0._EB
ALLOCATE(Q_DOT(N_Q_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','Q_DOT',IZERO)
Q_DOT = 0._EB
ALLOCATE(Q_DOT_SUM(N_Q_DOT,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','Q_DOT_SUM',IZERO)
Q_DOT_SUM = 0._EB
ALLOCATE(M_DOT(N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','M_DOT',IZERO)
M_DOT = 0._EB
ALLOCATE(M_DOT_SUM(N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','M_DOT_SUM',IZERO)
M_DOT_SUM=0._EB

ALLOCATE(MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','MASS_DT',IZERO)
MASS_DT=0._EB

ALLOCATE(PRESSURE_ERROR_MAX(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ERROR_MAX',IZERO)
ALLOCATE(PRESSURE_ERROR_MAX_LOC(3,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','PRESSURE_ERROR_MAX_LOC',IZERO)
PRESSURE_ERROR_MAX     = 0._EB
PRESSURE_ERROR_MAX_LOC = 0

ALLOCATE(VELOCITY_ERROR_MAX(NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX',IZERO)
ALLOCATE(VELOCITY_ERROR_MAX_LOC(3,NMESHES),STAT=IZERO)
CALL ChkMemErr('INIT','VELOCITY_ERROR_MAX_LOC',IZERO)
VELOCITY_ERROR_MAX     = 0._EB
VELOCITY_ERROR_MAX_LOC = 0

END SUBROUTINE INITIALIZE_GLOBAL_VARIABLES


!> \brief Initialize wall cell variables at external and obstruction boundaries
!> \param NM Mesh number
!> \param I x-index of inside wall cell
!> \param J y-index of inside wall cell
!> \param K z-index of inside wall cell
!> \param OBST_INDEX Index of the obstruction to which the wall cell is attached
!> \param IW Index of the wall cell
!> \param IOR Orientation index of the wall cell
!> \param SURF_INDEX Surface index of the wall cell
!> \param IERR Error code
!> \param TT Current time (s)

SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT
USE CONTROL_VARIABLES, ONLY : CONTROL
USE DEVICE_VARIABLES, ONLY : DEVICE
INTEGER, INTENT(IN) :: I,J,K,NM,OBST_INDEX,IW,IOR,SURF_INDEX
INTEGER  :: NOM_FOUND,NOM=0,ITER,IIO_MIN,IIO_MAX,JJO_MIN,JJO_MAX,KKO_MIN,KKO_MAX,VENT_INDEX
INTEGER, INTENT(OUT) :: IERR
REAL(EB), INTENT(IN) :: TT
REAL(EB) :: PX,PY,PZ,T_ACTIVATE,XIN,YIN,ZIN,DIST,XW,YW,ZW,RDN,AW,TSI,&
            ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,R1,RR,DELTA
INTEGER  :: N,SURF_INDEX_NEW,IIG,JJG,KKG,IIO,JJO,KKO,IC,ICG,ICO,NOM_CHECK(0:1),BOUNDARY_TYPE,FI
LOGICAL :: VENT_FOUND,ALIGNED
TYPE (MESH_TYPE), POINTER :: M,MM
TYPE (OBSTRUCTION_TYPE), POINTER :: OBX
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (WALL_TYPE), POINTER :: WC
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE (BOUNDARY_PROP2_TYPE), POINTER :: B2
TYPE (SURFACE_TYPE), POINTER :: SF

IERR = 0
M=>MESHES(NM)

! Determine if a VENT covers the surface

VENT_INDEX = 0
SURF_INDEX_NEW = SURF_INDEX
VENT_FOUND = .FALSE.

VENT_SEARCH_LOOP: DO N=1,M%N_VENT

   VT => M%VENTS(N)
   IF (OBST_INDEX>0) THEN
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY)           CYCLE VENT_SEARCH_LOOP
      IF (.NOT.M%OBSTRUCTION(OBST_INDEX)%ALLOW_VENT) CYCLE VENT_SEARCH_LOOP
      IF (VT%OBST_INDEX>0 .AND. VT%OBST_INDEX/=OBST_INDEX) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (VT%IOR/=IOR) CYCLE VENT_SEARCH_LOOP

   IF (ABS(IOR)==1) THEN
      IF (IOR== 1 .AND. I/=VT%I1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-1 .AND. I/=VT%I1+1) CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==2) THEN
      IF (IOR== 2 .AND. J/=VT%J1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-2 .AND. J/=VT%J1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (K<VT%K1+1 .OR. K>VT%K2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF
   IF (ABS(IOR)==3) THEN
      IF (IOR== 3 .AND. K/=VT%K1  ) CYCLE VENT_SEARCH_LOOP
      IF (IOR==-3 .AND. K/=VT%K1+1) CYCLE VENT_SEARCH_LOOP
      IF (I<VT%I1+1 .OR. I>VT%I2)   CYCLE VENT_SEARCH_LOOP
      IF (J<VT%J1+1 .OR. J>VT%J2)   CYCLE VENT_SEARCH_LOOP
      IF ( VT%RADIUS>0._EB .AND. ((M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2)>(VT%RADIUS**2) ) CYCLE VENT_SEARCH_LOOP
   ENDIF

   ! Check if there are over-lapping VENTs

   IF (VENT_FOUND) THEN
      WRITE(LU_ERR,'(A,I0,A,3(I0,1X),A,I0,A,I0,A)') 'WARNING: Two VENTs overlap in MESH ',NM,', Cell ',I,J,K,&
                                                    '. IOR ',IOR,'. VENT ',VT%ORDINAL,' rejected for that cell'
      EXIT VENT_SEARCH_LOOP
   ENDIF

   VENT_FOUND = .TRUE.

   ! Reassign the SURF index to be that of the VENT

   VENT_INDEX = N
   SURF_INDEX_NEW = VT%SURF_INDEX

ENDDO VENT_SEARCH_LOOP

! Compute boundary cell physical coords (XW,YW,ZW) and area (AW)

IF (ABS(IOR)==1) THEN
   IF (IOR== 1) THEN
      XW     = M%X(I)
      IIG    = I+1
      RDN    = M%RDXN(I)
      AW     = M%R(I)*M%DY(J)*M%DZ(K)
   ENDIF
   IF (IOR==-1) THEN
      XW     = M%X(I-1)
      IIG    = I-1
      RDN    = M%RDXN(I-1)
      AW     = M%R(I-1)*M%DY(J)*M%DZ(K)
   ENDIF
   JJG    = J
   KKG    = K
   YW = M%YC(J)
   ZW = M%ZC(K)
ENDIF
IF (ABS(IOR)==2) THEN
   IF (IOR== 2) THEN
      YW     = M%Y(J)
      JJG    = J+1
      RDN    = M%RDYN(J)
   ENDIF
   IF (IOR==-2) THEN
      YW     = M%Y(J-1)
      JJG    = J-1
      RDN    = M%RDYN(J-1)
   ENDIF
   IIG    = I
   KKG    = K
   XW = M%XC(I)
   ZW = M%ZC(K)
   AW = M%DX(I)*M%DZ(K)
ENDIF
IF (ABS(IOR)==3) THEN
   IF (IOR== 3) THEN
      ZW     = M%Z(K)
      KKG    = K+1
      RDN    = M%RDZN(K)
   ENDIF
   IF (IOR==-3) THEN
      ZW     = M%Z(K-1)
      KKG    = K-1
      RDN    = M%RDZN(K-1)
   ENDIF
   IIG    = I
   JJG    = J
   XW = M%XC(I)
   YW = M%YC(J)
   AW = M%DX(I)*M%RC(I)*M%DY(J)
ENDIF

IF (IOR==0) THEN
   IIG = I
   JJG = J
   KKG = K
ENDIF

! Save the wall index

IC  = M%CELL_INDEX(I  ,J  ,K  )
ICG = M%CELL_INDEX(IIG,JJG,KKG)

! Use BOUNDARY_TYPE to indicate whether the boundary cell is blocked or on an obstruction that is HIDDEN

BOUNDARY_TYPE = NULL_BOUNDARY

IF (IW<=M%N_EXTERNAL_WALL_CELLS .AND. OBST_INDEX==0) BOUNDARY_TYPE = SOLID_BOUNDARY

IF (OBST_INDEX>0) THEN
   IF (.NOT.M%OBSTRUCTION(OBST_INDEX)%HIDDEN) THEN
      BOUNDARY_TYPE = SOLID_BOUNDARY
      IF (IW<=M%N_EXTERNAL_WALL_CELLS) M%CELL(IC)%SOLID = .TRUE.
   ENDIF
ENDIF

IF (M%CELL(ICG)%SOLID) BOUNDARY_TYPE = NULL_BOUNDARY

! Check for neighboring meshes in a multiple mesh calculation

NOM_FOUND = 0
IIO_MIN =  1000000
IIO_MAX = -1000000
JJO_MIN =  1000000
JJO_MAX = -1000000
KKO_MIN =  1000000
KKO_MAX = -1000000
NOM_CHECK = 0

CHECK_MESHES: IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN

   DO ITER=0,1
      XIN = XW
      YIN = YW
      ZIN = ZW
      IF (SURF_INDEX_NEW==PERIODIC_SURF_INDEX .OR. SURF_INDEX_NEW==PERIODIC_FLOW_ONLY_SURF_INDEX) THEN
         SELECT CASE(IOR)
            CASE( 1) ; XIN = XF_MAX
            CASE(-1) ; XIN = XS_MIN
            CASE( 2) ; YIN = YF_MAX
            CASE(-2) ; YIN = YS_MIN
            CASE( 3) ; ZIN = ZF_MAX
            CASE(-3) ; ZIN = ZS_MIN
         END SELECT
      ENDIF
      IF (ABS(IOR)/=1) XIN = XW + (ITER*0.95_EB-0.475_EB)*(M%X(I)-M%X(I-1))
      IF (ABS(IOR)/=2) YIN = YW + (ITER*0.95_EB-0.475_EB)*(M%Y(J)-M%Y(J-1))
      IF (ABS(IOR)/=3) ZIN = ZW + (ITER*0.95_EB-0.475_EB)*(M%Z(K)-M%Z(K-1))
      IF (IOR== 1) XIN = XIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-1) XIN = XIN + MESH_SEPARATION_DISTANCE
      IF (IOR== 2) YIN = YIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-2) YIN = YIN + MESH_SEPARATION_DISTANCE
      IF (IOR== 3) ZIN = ZIN - MESH_SEPARATION_DISTANCE
      IF (IOR==-3) ZIN = ZIN + MESH_SEPARATION_DISTANCE
      CALL SEARCH_OTHER_MESHES(XIN,YIN,ZIN,NOM,IIO,JJO,KKO)
      NOM_CHECK(ITER) = NOM
      IF (NOM/=0) THEN
         IIO_MIN = MIN(IIO_MIN,IIO)
         IIO_MAX = MAX(IIO_MAX,IIO)
         JJO_MIN = MIN(JJO_MIN,JJO)
         JJO_MAX = MAX(JJO_MAX,JJO)
         KKO_MIN = MIN(KKO_MIN,KKO)
         KKO_MAX = MAX(KKO_MAX,KKO)
      ENDIF
   ENDDO

   ! Check to see if the current interpolated cell face spans more than one other mesh

   IF (NOM_CHECK(0)/=NOM_CHECK(1)) THEN
      WRITE(LU_ERR,'(A,I0,A,I0)') 'ERROR(431): MESH ',NM,' is not in alignment with MESH ',MAXVAL(NOM_CHECK)
      STOP_STATUS = SETUP_STOP
      IERR = 1
      RETURN
   ENDIF

   FOUND_OTHER_MESH: IF (NOM>0) THEN
      MM=>MESHES(NOM)
      ALIGNED = .TRUE.
      IF ( (ABS(IOR)==2 .OR. ABS(IOR)==3) .AND. MM%DX(IIO_MIN)<=M%DX(I) ) THEN
         IF (ABS( ((MM%X(IIO_MAX)-MM%X(IIO_MIN-1))-(M%X(I)-M%X(I-1))) / MM%DX(IIO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==3) .AND. MM%DY(JJO_MIN)<=M%DY(J) ) THEN
         IF (ABS( ((MM%Y(JJO_MAX)-MM%Y(JJO_MIN-1))-(M%Y(J)-M%Y(J-1))) / MM%DY(JJO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF ( (ABS(IOR)==1 .OR. ABS(IOR)==2) .AND. MM%DZ(KKO_MIN)<=M%DZ(K) ) THEN
         IF (ABS( ((MM%Z(KKO_MAX)-MM%Z(KKO_MIN-1))-(M%Z(K)-M%Z(K-1))) / MM%DZ(KKO_MIN))>ALIGNMENT_TOLERANCE ) ALIGNED = .FALSE.
      ENDIF
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I0,A,I0)') 'ERROR(431): MESH ',NM,' is out of alignment with MESH ',NOM
         STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      SELECT CASE(ABS(IOR))
         CASE(1)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(2)
            IF ( (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)>1.01_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) .AND. (M%DZ(K)<0.99_EB*MM%DZ(KKO_MIN)) ) ALIGNED = .FALSE.
         CASE(3)
            IF ( (M%DY(J)<0.99_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)>1.01_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
            IF ( (M%DY(J)>1.01_EB*MM%DY(JJO_MIN)) .AND. (M%DX(I)<0.99_EB*MM%DX(IIO_MIN)) ) ALIGNED = .FALSE.
      END SELECT
      IF (.NOT.ALIGNED) THEN
         WRITE(LU_ERR,'(A,I0,A,I0)') 'ERROR(431): MESH ',NM,' is out of alignment with MESH ',NOM
         STOP_STATUS = SETUP_STOP
         IERR = 1
         RETURN
      ENDIF

      ! NOM = "Number of the Other Mesh"
      !
      !     Mesh 1       Mesh 2
      !      3x6          1x2
      !  -------------------------
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   #->         |
      !  |---|---|---|           |
      !  |   |   | <-*           |
      !  |---|---|---|-----------|
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   |           |
      !  |---|---|---|           |
      !  |   |   |   |           |
      !  -------------------------
      !
      !                         NM=1,NOM=2,IW=*        NM=2,NOM=1,IW=#

      NOM_FOUND = NOM
      ICO = MM%CELL_INDEX(IIO_MIN,JJO_MIN,KKO_MIN)

      IF (OBST_INDEX==0) THEN
         IF (.NOT.M%CELL(ICG)%SOLID .AND. .NOT.MM%CELL(ICO)%SOLID) THEN
            BOUNDARY_TYPE = INTERPOLATED_BOUNDARY
            IF (SURF_INDEX_NEW/=PERIODIC_FLOW_ONLY_SURF_INDEX) SURF_INDEX_NEW = INTERPOLATED_SURF_INDEX
         ENDIF
         IF (M%CELL(ICG)%SOLID .OR. MM%CELL(ICO)%SOLID) THEN
            IF (MM%CELL(ICO)%SOLID) THEN
               SURF_INDEX_NEW = MM%OBSTRUCTION(MM%CELL(ICO)%OBST_INDEX)%SURF_INDEX(IOR)
            ELSE
               SURF_INDEX_NEW = M%OBSTRUCTION(M%CELL(ICG)%OBST_INDEX)%SURF_INDEX(-IOR)
            ENDIF
         ENDIF
         VENT_INDEX = 0
      ENDIF

      ! Determine if masses from consumable obstructions need to be exchanged

      IF (M%CELL(ICG)%SOLID .OR. MM%CELL(ICO)%SOLID) THEN
         IF (M%OBSTRUCTION(M%CELL(ICG)%OBST_INDEX)%CONSUMABLE .OR. MM%OBSTRUCTION(MM%CELL(ICO)%OBST_INDEX)%CONSUMABLE) &
            EXCHANGE_OBST_MASS = .TRUE.
      ENDIF

      ! Do not allow a MIRROR boundary to sit on a mesh interface

      IF (VENT_INDEX>0) THEN
         IF (M%VENTS(VENT_INDEX)%BOUNDARY_TYPE==MIRROR_BOUNDARY) VENT_INDEX = 0
      ENDIF

      ! Open up the ghost cell at the interpolated boundary

      IF (BOUNDARY_TYPE == INTERPOLATED_BOUNDARY) M%CELL(M%CELL_INDEX(I,J,K))%SOLID = .FALSE.

   ENDIF FOUND_OTHER_MESH

ENDIF CHECK_MESHES

! If this wall cell is attached to an OBST, check if the OBST face is exposed

IF (OBST_INDEX>0) THEN
   IF (.NOT.M%CELL(ICG)%SOLID .OR. M%OBSTRUCTION(M%CELL(ICG)%OBST_INDEX)%REMOVABLE) THEN
      FI = ABS(IOR)*2 ; IF (IOR<0) FI = FI-1
      M%OBSTRUCTION(OBST_INDEX)%EXPOSED_FACE_INDEX(FI) = 1
   ENDIF
ENDIF

! Ensure that the WALL_INDEX and SURF_INDEX can be identified from the abutting gas phase cell, ICG

M%CELL(ICG)%WALL_INDEX(-IOR) = IW
M%CELL(ICG)%SURF_INDEX(-IOR) = SURF_INDEX_NEW

! Ensure that there is an open slot in M%WALL and its associated derived types

CALL ALLOCATE_STORAGE(NM,WALL_INDEX=IW,SURF_INDEX=SURF_INDEX_NEW)

! Initialize wall cell (WC) variables

SF => SURFACE(SURF_INDEX_NEW)
WC => M%WALL(IW)

WC%SURF_INDEX = SURF_INDEX_NEW
WC%OBST_INDEX = OBST_INDEX
WC%BOUNDARY_TYPE = BOUNDARY_TYPE

IF (IW<=M%N_EXTERNAL_WALL_CELLS) THEN
   EWC => M%EXTERNAL_WALL(IW)
   EWC%NOM     = NOM_FOUND
   EWC%IIO_MIN = IIO_MIN
   EWC%JJO_MIN = JJO_MIN
   EWC%KKO_MIN = KKO_MIN
   EWC%IIO_MAX = IIO_MAX
   EWC%JJO_MAX = JJO_MAX
   EWC%KKO_MAX = KKO_MAX
ENDIF

BC => M%BOUNDARY_COORD(WC%BC_INDEX)

BC%II = I
BC%JJ = J
BC%KK = K
BC%IIG = IIG
BC%JJG = JJG
BC%KKG = KKG
BC%IOR = IOR
SELECT CASE(BC%IOR)
   CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/)
   CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/)
   CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/)
   CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/)
   CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/)
   CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/)
END SELECT
BC%X = XW
BC%Y = YW
BC%Z = ZW
SELECT CASE(BC%IOR)
   CASE(-1) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I-1) ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
   CASE( 1) ; BC%X1=M%X(I)   ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
   CASE(-2) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J-1) ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
   CASE( 2) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J)   ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
   CASE(-3) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K-1)
   CASE( 3) ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)   ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)   ; BC%Z1=M%Z(K)   ; BC%Z2=M%Z(K)
END SELECT

! If the WALL cell is attached to a THIN obstruction, use the obstruction coords for the wall cell coordinates

IF (OBST_INDEX>0) THEN
   OBX => M%OBSTRUCTION(OBST_INDEX)
   IF (OBX%I1==OBX%I2 .AND. OBX%THIN .AND. .NOT.M%CELL(IC)%SOLID) THEN ; BC%X1=OBX%X1 ; BC%X2=OBX%X2 ; ENDIF
   IF (OBX%J1==OBX%J2 .AND. OBX%THIN .AND. .NOT.M%CELL(IC)%SOLID) THEN ; BC%Y1=OBX%Y1 ; BC%Y2=OBX%Y2 ; ENDIF
   IF (OBX%K1==OBX%K2 .AND. OBX%THIN .AND. .NOT.M%CELL(IC)%SOLID) THEN ; BC%Z1=OBX%Z1 ; BC%Z2=OBX%Z2 ; ENDIF
ENDIF

B1 => M%BOUNDARY_PROP1(WC%B1_INDEX)
B2 => M%BOUNDARY_PROP2(WC%B2_INDEX)

B2%U_TAU = 0._EB
B2%Y_PLUS = 1._EB
B2%Z_STAR = 1._EB
B2%HEAT_TRANSFER_REGIME = 0

B1%RDN    = RDN
B1%AREA   = AW

! If the simulation is only a TGA analysis, get the wall index

IF (WC%SURF_INDEX==TGA_SURF_INDEX) THEN
   TGA_WALL_INDEX = IW
   TGA_MESH_INDEX = NM
ENDIF

! Assign internal values of temp, density, and mass fraction

B1%RHO_F = M%RHO(IIG,JJG,KKG)
B1%RHO_D_F = 0._EB
B1%RHO_D_DZDN_F = 0._EB

IF (M%CELL(M%CELL_INDEX(I,J,K))%SOLID .OR. IW<=M%N_EXTERNAL_WALL_CELLS) THEN
   M%RSUM(I,J,K) = M%RSUM(IIG,JJG,KKG)
   B1%ZZ_F(1:N_TRACKED_SPECIES)  = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
   M%ZZ(I,J,K,1:N_TRACKED_SPECIES) = M%ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES)
ENDIF

! Compute the mass of the grid cell corresponding to the wall cell

IF (OBST_INDEX>0) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   IF (OBX%FDS_AREA(ABS(IOR))>TWO_EPSILON_EB) B1%AREA_ADJUST = SF%AREA_MULTIPLIER*OBX%INPUT_AREA(ABS(IOR))/OBX%FDS_AREA(ABS(IOR))
   IF (B1%AREA_ADJUST<=TWO_EPSILON_EB) B1%AREA_ADJUST = 1._EB
   IF (OBX%MASS>1.E5_EB .AND. OBX%BULK_DENSITY<=0._EB) OBX%MASS = SF%SURFACE_DENSITY*B1%AREA*B1%AREA_ADJUST
ENDIF

! Prescribe normal velocity for surface cell (U_NORMAL_0)

B1%U_NORMAL_0 = SF%VEL

IF (OBST_INDEX>0 .AND. ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   IF (OBX%FDS_AREA(ABS(IOR))>TWO_EPSILON_EB) &
      B1%U_NORMAL_0 = SF%VOLUME_FLOW*(OBX%INPUT_AREA(ABS(IOR))/OBX%UNDIVIDED_INPUT_AREA(ABS(IOR))) / OBX%FDS_AREA(ABS(IOR))
ENDIF
IF (OBST_INDEX>0 .AND. ABS(SF%MASS_FLUX_TOTAL)>=TWO_EPSILON_EB) THEN
   OBX=>M%OBSTRUCTION(OBST_INDEX)
   B1%U_NORMAL_0 = SF%MASS_FLUX_TOTAL / RHOA * B1%AREA_ADJUST
ENDIF

! Do VENT-specific set-ups

T_ACTIVATE    = T_BEGIN-1._EB
WC%VENT_INDEX = VENT_INDEX

PROCESS_VENT: IF (WC%VENT_INDEX>0) THEN

   VT => M%VENTS(WC%VENT_INDEX)

   B1%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%INPUT_AREA/VT%FDS_AREA
   IF (B1%AREA_ADJUST<=TWO_EPSILON_EB) B1%AREA_ADJUST = 1._EB

   IF (VT%CTRL_INDEX > 0) THEN
      IF (.NOT. CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) T_ACTIVATE=1E10_EB
   ENDIF
   IF (VT%DEVC_INDEX > 0) THEN
      IF (.NOT. DEVICE(VT%DEVC_INDEX)%CURRENT_STATE) T_ACTIVATE=1E10_EB
   ENDIF

   ! Set the velocity at each surface cell

   B1%U_NORMAL_0 = SF%VEL

   IF (ABS(SF%VOLUME_FLOW)>TWO_EPSILON_EB)   THEN
      B1%U_NORMAL_0 = SF%VOLUME_FLOW*(VT%INPUT_AREA/VT%UNDIVIDED_INPUT_AREA)/VT%FDS_AREA
   ENDIF
   IF (ABS(SF%MASS_FLUX_TOTAL)>TWO_EPSILON_EB) B1%U_NORMAL_0 = SF%MASS_FLUX_TOTAL/RHOA*B1%AREA_ADJUST

   IF (SF%CONVERT_VOLUME_TO_MASS) THEN
      IF (ABS(B1%U_NORMAL_0)>TWO_EPSILON_EB) THEN
         ZZ_GET=0._EB
         ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,SF%MASS_FRACTION(1:N_TRACKED_SPECIES))
         CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
         SF%MASS_FLUX = -RHOA*(RSUM0/RSUM_F)*(TMPA/SF%TMP_FRONT)*SF%MASS_FRACTION*B1%U_NORMAL_0
         SF%SPECIES_BC_INDEX = SPECIFIED_MASS_FLUX
      ELSE
         CALL SHUTDOWN('ERROR(432): SURF: '//TRIM(SF%ID)//' must specify velocity boundary condition for conversion',&
                        PROCESS_0_ONLY=.FALSE.)
         IERR = 1
         RETURN
      ENDIF
   ENDIF

   ! Special velocity profiles

   PARABOLIC_IF: IF (SF%PROFILE==PARABOLIC_PROFILE) THEN
      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2
               B1%U_NORMAL_0 = B1%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PY = 4._EB*(M%YC(J)-VT%Y1_ORIG)*(VT%Y2_ORIG-M%YC(J))/(VT%Y2_ORIG-VT%Y1_ORIG)**2
               PZ = 4._EB*(M%ZC(K)-VT%Z1_ORIG)*(VT%Z2_ORIG-M%ZC(K))/(VT%Z2_ORIG-VT%Z1_ORIG)**2
               B1%U_NORMAL_0 = B1%U_NORMAL_0*PY*PZ
            ENDIF
         CASE(2)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2
               B1%U_NORMAL_0 = B1%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PX = 4._EB*(M%XC(I)-VT%X1_ORIG)*(VT%X2_ORIG-M%XC(I))/(VT%X2_ORIG-VT%X1_ORIG)**2
               PZ = 4._EB*(M%ZC(K)-VT%Z1_ORIG)*(VT%Z2_ORIG-M%ZC(K))/(VT%Z2_ORIG-VT%Z1_ORIG)**2
               B1%U_NORMAL_0 = B1%U_NORMAL_0*PX*PZ
            ENDIF
         CASE(3)
            IF (VT%RADIUS>0._EB) THEN
               RR = (M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2
               B1%U_NORMAL_0 = B1%U_NORMAL_0*(VT%RADIUS**2-RR)/VT%RADIUS**2
            ELSE
               PX = 4._EB*(M%XC(I)-VT%X1_ORIG)*(VT%X2_ORIG-M%XC(I))/(VT%X2_ORIG-VT%X1_ORIG)**2
               PY = 4._EB*(M%YC(J)-VT%Y1_ORIG)*(VT%Y2_ORIG-M%YC(J))/(VT%Y2_ORIG-VT%Y1_ORIG)**2
               IF (CYLINDRICAL) THEN
                  B1%U_NORMAL_0 = B1%U_NORMAL_0*PX
               ELSE
                  B1%U_NORMAL_0 = B1%U_NORMAL_0*PX*PY
               ENDIF
            ENDIF
      END SELECT
      IF (ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB) THEN   ! Match desired volume flow
         IF (VT%RADIUS>0._EB) THEN
            B1%U_NORMAL_0 = B1%U_NORMAL_0*2._EB
         ELSE
            B1%U_NORMAL_0 = B1%U_NORMAL_0*9._EB/4._EB
         ENDIF
      ENDIF
   ENDIF PARABOLIC_IF

   IF (SF%PROFILE==BOUNDARY_LAYER_PROFILE) THEN

      ! Currently only set up for circular vents

      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%YC(J)-VT%Y0)**2 + (M%ZC(K)-VT%Z0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  B1%U_NORMAL_0 = B1%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
         CASE(2)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%XC(I)-VT%X0)**2 + (M%ZC(K)-VT%Z0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  B1%U_NORMAL_0 = B1%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
         CASE(3)
            IF (VT%RADIUS>0._EB) THEN
               DELTA = VT%RADIUS - SQRT( VT%RADIUS**2*(2._EB*ABS(SF%VEL_BULK/SF%VEL)-1._EB) )
               R1 = VT%RADIUS - DELTA
               RR = SQRT( (M%XC(I)-VT%X0)**2 + (M%YC(J)-VT%Y0)**2 )
               IF (RR>R1 .AND. RR<=VT%RADIUS .AND. DELTA>TWO_EPSILON_EB) THEN
                  B1%U_NORMAL_0 = B1%U_NORMAL_0*(1._EB - ((RR-R1)/DELTA)**2 )
               ENDIF
            ENDIF
      END SELECT
   ENDIF

   IF (SF%PROFILE==ATMOSPHERIC_PROFILE) THEN
      IF (M%ZC(K)<GROUND_LEVEL) THEN
         CALL SHUTDOWN('ERROR(433): SURF '//TRIM(SF%ID)//' cannot be applied below GROUND_LEVEL.',PROCESS_0_ONLY=.FALSE.)
         IERR = 1
         RETURN
      ENDIF
      B1%U_NORMAL_0 = B1%U_NORMAL_0*((M%ZC(K)-GROUND_LEVEL)/SF%Z0)**SF%PLE
   ENDIF

   IF (SF%PROFILE==RAMP_PROFILE) THEN
      SELECT CASE(ABS(IOR))
         CASE(1)
            IF (SF%RAMP(VELO_PROF_X)%ID/='null') THEN
               CALL SHUTDOWN('ERROR(434): RAMP_V_X assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%YC(J),SF%RAMP(VELO_PROF_Y)%INDEX,TAU=1._EB)
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%ZC(K),SF%RAMP(VELO_PROF_Z)%INDEX,TAU=1._EB)
         CASE(2)
            IF (SF%RAMP(VELO_PROF_Y)%ID/='null') THEN
               CALL SHUTDOWN('ERROR(435): RAMP_V_Y assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%XC(I),SF%RAMP(VELO_PROF_X)%INDEX,TAU=1._EB)
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%ZC(K),SF%RAMP(VELO_PROF_Z)%INDEX,TAU=1._EB)
         CASE(3)
            IF (SF%RAMP(VELO_PROF_Z)%ID/='null') THEN
               CALL SHUTDOWN('ERROR(436): RAMP_V_Z assigned to SURF '//TRIM(SF%ID),PROCESS_0_ONLY=.FALSE.)
               IERR = 1
               RETURN
            ENDIF
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%XC(I),SF%RAMP(VELO_PROF_X)%INDEX,TAU=1._EB)
            B1%U_NORMAL_0 = B1%U_NORMAL_0*EVALUATE_RAMP(M%YC(J),SF%RAMP(VELO_PROF_Y)%INDEX,TAU=1._EB)
      END SELECT
   ENDIF

   ! Check if fire spreads radially over this vent

   IF (VT%FIRE_SPREAD_RATE>0._EB) THEN
      DIST = SQRT((BC%X-VT%X0)**2 + (BC%Y-VT%Y0)**2 + (BC%Z-VT%Z0)**2)
      T_ACTIVATE = TT + DIST/VT%FIRE_SPREAD_RATE
   ENDIF

   ! Miscellaneous settings

   IF (.NOT.M%CELL(ICG)%SOLID) THEN
      IF (VT%BOUNDARY_TYPE==MIRROR_BOUNDARY) THEN
         WC%BOUNDARY_TYPE = MIRROR_BOUNDARY
         WC%SURF_INDEX    = MIRROR_SURF_INDEX
      ENDIF
      IF (VT%BOUNDARY_TYPE==OPEN_BOUNDARY) THEN
         WC%BOUNDARY_TYPE = OPEN_BOUNDARY
         WC%SURF_INDEX    = OPEN_SURF_INDEX
      ENDIF
   ENDIF

ENDIF PROCESS_VENT

! Check if fire spreads radially over this surface type

IF (SF%FIRE_SPREAD_RATE>0._EB) THEN
   DIST = SQRT((BC%X-SF%XYZ(1))**2 +(BC%Y-SF%XYZ(2))**2 +(BC%Z-SF%XYZ(3))**2)
   T_ACTIVATE = TT + DIST/SF%FIRE_SPREAD_RATE
ENDIF

! Set ignition time of each boundary cell

IF (T_ACTIVATE < T_BEGIN) THEN
   IF (SF%T_IGN==T_BEGIN) THEN
      B1%T_IGN = TT
   ELSE
      B1%T_IGN = SF%T_IGN
   ENDIF
ELSE
   B1%T_IGN = T_ACTIVATE
ENDIF

! Set correct initial value of temperature for RAMP_T

IF (ABS(B1%T_IGN-T_BEGIN) <= SPACING(B1%T_IGN) .AND. SF%RAMP(TIME_TEMP)%INDEX>=1) THEN
   TSI = TT
ELSE
   TSI = TT - B1%T_IGN
ENDIF

IF (SF%RAMP_T_I_INDEX > 0) THEN
   B1%TMP_F = EVALUATE_RAMP(0._EB,SF%RAMP_T_I_INDEX)
   B1%TMP_B = EVALUATE_RAMP(SUM(SF%LAYER_THICKNESS),SF%RAMP_T_I_INDEX)
ELSE
   IF (SF%TMP_FRONT_INITIAL>0._EB) THEN
      B1%TMP_F = SF%TMP_FRONT_INITIAL
   ELSEIF (SF%TMP_FRONT>0._EB) THEN
      B1%TMP_F = M%TMP_0(BC%KK) + &
                 EVALUATE_RAMP(TSI,SF%RAMP(TIME_TEMP)%INDEX,TAU=SF%RAMP(TIME_TEMP)%TAU)*(SF%TMP_FRONT-M%TMP_0(BC%KK))
   ELSE
      B1%TMP_F = M%TMP_0(BC%KK)
   ENDIF

   IF (SF%TMP_BACK>0._EB) THEN
      B1%TMP_B = SF%TMP_BACK
   ELSE
      B1%TMP_B = SF%TMP_INNER
   ENDIF
ENDIF

! Reinitialize wall cell outgoing radiation for change in TMP_F

IF (RADIATION) B1%Q_RAD_OUT = SF%EMISSIVITY*SIGMA*B1%TMP_F**4

! Record original boundary condition index for exterior wall cells that might get covered up

IF (OBST_INDEX==0 .AND. IW<=M%N_EXTERNAL_WALL_CELLS) EWC%SURF_INDEX_ORIG = SURF_INDEX_NEW

END SUBROUTINE INIT_WALL_CELL


!> \brief Initialize thin wall cell variables at edges of thin obstructions when 3-D heat transfer is specified
!> \param NM Mesh number
!> \param I x-index of inside thin wall cell
!> \param J y-index of inside thin wall cell
!> \param K z-index of inside thin wall cell
!> \param OBST_INDEX Index of the obstruction to which the thin wall cell is attached
!> \param ITW Index of the thin wall cell
!> \param IOR Orientation index of the thin wall cell
!> \param SURF_INDEX Surface index of the thin wall cell
!> \param IEC Edge index

SUBROUTINE INIT_THIN_WALL_CELL(NM,I,J,K,OBST_INDEX,ITW,IOR,SURF_INDEX,IEC)

USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN) :: I,J,K,NM,OBST_INDEX,ITW,IOR,SURF_INDEX,IEC
INTEGER :: IC
TYPE (MESH_TYPE), POINTER :: M
TYPE (THIN_WALL_TYPE), POINTER :: TW
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (OBSTRUCTION_TYPE), POINTER :: OB

M=>MESHES(NM)

OB => M%OBSTRUCTION(OBST_INDEX)
IC = M%CELL_INDEX(I,J,K)
M%CELL(IC)%THIN_WALL_INDEX(IOR,IEC) = ITW
M%CELL(IC)%THIN_SURF_INDEX(IOR,IEC) = SURF_INDEX
M%CELL(IC)%THIN_OBST_INDEX(IOR,IEC) = OB%ORDINAL

! Ensure that there is an open slot in M%WALL and its associated derived types

CALL ALLOCATE_STORAGE(NM,THIN_WALL_INDEX=ITW,SURF_INDEX=SURF_INDEX)

! Initialize thin wall cell (TW) variables

SF => SURFACE(SURF_INDEX)
TW => M%THIN_WALL(ITW)

TW%SURF_INDEX = SURF_INDEX
TW%OBST_INDEX = OBST_INDEX
TW%BOUNDARY_TYPE = SOLID_BOUNDARY
TW%IEC = IEC

BC => M%BOUNDARY_COORD(TW%BC_INDEX)

BC%II = I
BC%JJ = J
BC%KK = K
BC%IIG = I
BC%JJG = J
BC%KKG = K
BC%IOR = IOR
SELECT CASE(BC%IOR)
   CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/)
   CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/)
   CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/)
   CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/)
   CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/)
   CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/)
END SELECT
BC%X  = M%X(I)
BC%Y  = M%Y(J)
BC%Z  = M%Z(K)
SELECT CASE(ABS(BC%IOR))
   CASE(1)
      SELECT CASE(BC%IOR)
         CASE(-1) ; BC%X1=OB%X1 ; BC%X2=OB%X1
         CASE( 1) ; BC%X1=OB%X2 ; BC%X2=OB%X2
      END SELECT
      SELECT CASE(IEC)
         CASE(2) ; BC%Z1=OB%Z1 ; BC%Z2=OB%Z2 ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)
         CASE(3) ; BC%Y1=OB%Y1 ; BC%Y2=OB%Y2 ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
      END SELECT
   CASE(2)
      SELECT CASE(BC%IOR)
         CASE(-2) ; BC%Y1=OB%Y1 ; BC%Y2=OB%Y1
         CASE( 2) ; BC%Y1=OB%Y2 ; BC%Y2=OB%Y2
      END SELECT
      SELECT CASE(IEC)
         CASE(1) ; BC%Z1=OB%Z1 ; BC%Z2=OB%Z2 ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)
         CASE(3) ; BC%X1=OB%X1 ; BC%X2=OB%X2 ; BC%Z1=M%Z(K-1) ; BC%Z2=M%Z(K)
      END SELECT
   CASE(3)
      SELECT CASE(BC%IOR)
         CASE(-3) ; BC%Z1=OB%Z1 ; BC%Z2=OB%Z1
         CASE( 3) ; BC%Z1=OB%Z2 ; BC%Z2=OB%Z2
      END SELECT
      SELECT CASE(IEC)
         CASE(1) ; BC%Y1=OB%Y1 ; BC%Y2=OB%Y2 ; BC%X1=M%X(I-1) ; BC%X2=M%X(I)
         CASE(2) ; BC%X1=OB%X1 ; BC%X2=OB%X2 ; BC%Y1=M%Y(J-1) ; BC%Y2=M%Y(J)
      END SELECT
END SELECT

B1 => M%BOUNDARY_PROP1(TW%B1_INDEX)

IF (SF%RAMP_T_I_INDEX > 0) THEN
   B1%TMP_F = EVALUATE_RAMP(0._EB,SF%RAMP_T_I_INDEX)
   B1%TMP_B = EVALUATE_RAMP(SUM(SF%LAYER_THICKNESS),SF%RAMP_T_I_INDEX)
ELSE
   IF (SF%TMP_FRONT_INITIAL>0._EB) THEN
      B1%TMP_F = SF%TMP_FRONT_INITIAL
   ELSEIF (SF%TMP_FRONT>0._EB) THEN
      B1%TMP_F = M%TMP_0(BC%KK) + &
                 EVALUATE_RAMP(T_BEGIN,SF%RAMP(TIME_TEMP)%INDEX,TAU=SF%RAMP(TIME_TEMP)%TAU)*(SF%TMP_FRONT-M%TMP_0(BC%KK))
   ELSE
      B1%TMP_F = M%TMP_0(BC%KK)
   ENDIF

   IF (SF%TMP_BACK>0._EB) THEN
      B1%TMP_B = SF%TMP_BACK
   ELSE
      B1%TMP_B = SF%TMP_INNER
   ENDIF
ENDIF

! Reinitialize wall cell outgoing radiation for change in TMP_F

IF (RADIATION) B1%Q_RAD_OUT = SF%EMISSIVITY*SIGMA*B1%TMP_F**4

END SUBROUTINE INIT_THIN_WALL_CELL


!> \brief Locate wall back indices
!> \param NM Mesh number
!> \details Loop through all internal and external wall cells and look for thermally-thick
!> solids with EXPOSED back wall cells. If the exposed back wall cell is in
!> another mesh, store the cell info into arrays that are to be MPI exchanged.

SUBROUTINE FIND_WALL_BACK_INDICES(NM)

USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_INTEGER_ARRAY
INTEGER, INTENT(IN) :: NM
INTEGER :: IW,ITW,N,NOM,IC,IOR,IEC
TYPE(MESH_TYPE), POINTER :: M,M4
TYPE(STORAGE_TYPE), POINTER :: OS

M => MESHES(NM)

! Find and save the back mesh and indices for all WALL cells in the current mesh

DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
   CALL FIND_WALL_BACK_INDEX(NM,IW)
ENDDO

! Search all neighboring meshes for 3-D WALL cells. Add index and surface information from these to M%OMESH(NOM)%WALL_RECV_BUFFER

DO N=1,M%N_NEIGHBORING_MESHES
   NOM = M%NEIGHBORING_MESH(N)
   IF (NM==NOM) CYCLE
   M4 => MESHES(NOM)
   IF ((M%XS>=M4%XF .OR. M%XF<=M4%XS) .AND. (M%YS>=M4%YF .OR. M%YF<=M4%YS) .AND. (M%ZS>=M4%ZF .OR. M%ZF<=M4%ZS)) CYCLE
   OS => M%OMESH(NOM)%WALL_RECV_BUFFER
   DO IC=1,CELL_COUNT(NOM)
      IF (M4%CELL(IC)%SOLID) CYCLE
      DO IOR=-3,3
         IF (IOR==0) CYCLE
         IF (SURFACE(M4%CELL(IC)%SURF_INDEX(IOR))%HT_DIM==1) CYCLE
         IF (.NOT.ALLOCATED(OS%ITEM_INDEX)) THEN
            OS%N_ITEMS_DIM = 50
            ALLOCATE(OS%ITEM_INDEX(1:OS%N_ITEMS_DIM))
            ALLOCATE(OS%SURF_INDEX(1:OS%N_ITEMS_DIM))
         ENDIF
         IF (COUNT(OS%ITEM_INDEX(1:OS%N_ITEMS)==M4%CELL(IC)%WALL_INDEX(IOR))>0) CYCLE
         IF (OS%N_ITEMS>=OS%N_ITEMS_DIM) THEN
            CALL REALLOCATE_INTEGER_ARRAY(OS%ITEM_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
            CALL REALLOCATE_INTEGER_ARRAY(OS%SURF_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
            OS%N_ITEMS_DIM = OS%N_ITEMS_DIM + 50
         ENDIF
         OS%N_ITEMS = OS%N_ITEMS + 1
         OS%ITEM_INDEX(OS%N_ITEMS) = M4%CELL(IC)%WALL_INDEX(IOR)
         OS%SURF_INDEX(OS%N_ITEMS) = M4%CELL(IC)%SURF_INDEX(IOR)
      ENDDO
   ENDDO
ENDDO

! Find back index of thin wall

DO ITW=1,M%N_THIN_WALL_CELLS
   CALL FIND_THIN_WALL_BACK_INDEX(NM,ITW)
ENDDO

! Search all neighboring meshes for 3-D THIN_WALL cells. Add index and surface info from these to M%OMESH(NOM)%THIN_WALL_RECV_BUFFER

DO N=1,M%N_NEIGHBORING_MESHES
   NOM = M%NEIGHBORING_MESH(N)
   IF (NM==NOM) CYCLE
   M4 => MESHES(NOM)
   OS => M%OMESH(NOM)%THIN_WALL_RECV_BUFFER
   DO IC=1,CELL_COUNT(NOM)
      DO IOR=-3,3
         IF (IOR==0) CYCLE
         DO IEC=1,3
            IF (M4%CELL(IC)%THIN_WALL_INDEX(IOR,IEC)>0) THEN
               IF (OS%N_ITEMS>0) THEN
                  IF (COUNT(OS%ITEM_INDEX(1:OS%N_ITEMS)==M4%CELL(IC)%THIN_WALL_INDEX(IOR,IEC))>0) CYCLE
               ENDIF
               IF (OS%N_ITEMS>=OS%N_ITEMS_DIM) THEN
                  CALL REALLOCATE_INTEGER_ARRAY(OS%ITEM_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
                  CALL REALLOCATE_INTEGER_ARRAY(OS%SURF_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
                  OS%N_ITEMS_DIM = OS%N_ITEMS_DIM + 50
               ENDIF
               OS%N_ITEMS = OS%N_ITEMS + 1
               OS%ITEM_INDEX(OS%N_ITEMS) = M4%CELL(IC)%THIN_WALL_INDEX(IOR,IEC)
               OS%SURF_INDEX(OS%N_ITEMS) = M4%CELL(IC)%THIN_SURF_INDEX(IOR,IEC)
            ENDIF
         ENDDO
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE FIND_WALL_BACK_INDICES


!> \brief Find the back wall cell for a given wall cell.
!> \param NM Mesh number
!> \param IW Wall cell index
!> \details If the exposed back wall cell is in another mesh, store the cell info into arrays that are to be MPI exchanged.

SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)

USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_INTEGER_ARRAY
USE MISC_FUNCTIONS, ONLY: PROCESS_MESH_NEIGHBORHOOD
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
INTEGER, INTENT(IN) :: NM,IW
INTEGER :: II,JJ,KK,IC,ICG,IOR,NOM,ITER,OBST_INDEX,OBST_INDEX_PREVIOUS,NN,NNN,NL,N_LAYERS_OBST,&
           N_MATL_OBST,N_LAYERS,N_MATLS,IIF,JJF,KKF
INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX_OBST,MATL_INDEX
REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION
REAL(EB), DIMENSION(0:MAX_LAYERS) :: LAYER_THICKNESS,MINIMUM_LAYER_THICKNESS
REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,&
                                   STRETCH_FACTOR,STRETCH_FACTOR_OBST,CELL_SIZE,CELL_SIZE_OBST,&
                                   CELL_SIZE_FACTOR,CELL_SIZE_FACTOR_OBST
INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_MAX,N_LAYER_CELLS_MAX_OBST,RAMP_IHS_INDEX,RAMP_IHS_INDEX_OBST
LOGICAL, DIMENSION(MAX_LAYERS) :: HT3D_LAYER
REAL(EB) :: XXC,YYC,ZZC,THICKNESS,OLD_THICKNESS,FRONT_LINING_THICKNESS,BACK_LINING_THICKNESS,LAYER_THICKNESS_OBST_TOTAL
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
LOGICAL :: THIN_OBSTRUCTION
TYPE (MESH_TYPE), POINTER :: M
TYPE (WALL_TYPE), POINTER :: WC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE (MESH_TYPE), POINTER :: OM,OM_PREV
TYPE (OBSTRUCTION_TYPE), POINTER :: OB,OB_PREV
TYPE (SURFACE_TYPE), POINTER :: SF,SF_BACK
TYPE (STORAGE_TYPE), POINTER :: OS
LOGICAL :: SUCCESS

M => MESHES(NM)
WC => M%WALL(IW)
SF => SURFACE(WC%SURF_INDEX)
IF (SF%THERMAL_BC_INDEX/=THERMALLY_THICK) RETURN
BC => M%BOUNDARY_COORD(WC%BC_INDEX)
IF (TWO_D .AND. (ABS(BC%IOR)==2.OR.(CYLINDRICAL.AND.BC%IOR==1)) .AND. IW<=M%N_EXTERNAL_WALL_CELLS) RETURN
ICG = M%CELL_INDEX(BC%IIG,BC%JJG,BC%KKG)
IF (M%CELL(ICG)%SOLID) RETURN
ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX)
II  = BC%II
JJ  = BC%JJ
KK  = BC%KK
IOR = BC%IOR
NOM = NM
OM => MESHES(NOM)
OM_PREV => MESHES(NOM)
ITER = 0
OBST_INDEX = 0
THICKNESS  = 0._EB
THIN_OBSTRUCTION = .FALSE.
MATL_INDEX_OBST = 0
MATL_INDEX = 0
IF (SF%VARIABLE_THICKNESS .OR. SF%HT_DIM>1) THEN
   N_LAYERS_OBST = 1
   LAYER_THICKNESS_OBST = 0._EB
   MATL_MASS_FRACTION_OBST = 0._EB
   N_MATL_OBST = 0
   HEAT_SOURCE_OBST = 0._EB
   RAMP_IHS_INDEX_OBST = -1
   STRETCH_FACTOR_OBST = SF%STRETCH_FACTOR(1)
   CELL_SIZE_OBST = SF%CELL_SIZE(1)
   CELL_SIZE_FACTOR_OBST = SF%CELL_SIZE_FACTOR(1)
   N_LAYER_CELLS_MAX_OBST = SF%N_LAYER_CELLS_MAX(1)
ENDIF

FIND_BACK_WALL_CELL: DO  ! Look for the back wall cell; that is, the wall cell on the other side of the obstruction

   ITER = ITER + 1
   OM_PREV => MESHES(NOM)

   IF (II==0 .OR. II==OM%IBP1 .OR. JJ==0 .OR. JJ==OM%JBP1 .OR. KK==0 .OR. KK==OM%KBP1) THEN
      XXC=OM%XC(II) ; YYC=OM%YC(JJ) ; ZZC=OM%ZC(KK)
      IF (II==0)       XXC = OM%X(II)   - MESH_SEPARATION_DISTANCE
      IF (II==OM%IBP1) XXC = OM%X(II-1) + MESH_SEPARATION_DISTANCE
      IF (JJ==0)       YYC = OM%Y(JJ)   - MESH_SEPARATION_DISTANCE
      IF (JJ==OM%JBP1) YYC = OM%Y(JJ-1) + MESH_SEPARATION_DISTANCE
      IF (KK==0)       ZZC = OM%Z(KK)   - MESH_SEPARATION_DISTANCE
      IF (KK==OM%KBP1) ZZC = OM%Z(KK-1) + MESH_SEPARATION_DISTANCE
      CALL SEARCH_OTHER_MESHES(XXC,YYC,ZZC,NOM,II,JJ,KK)
      IF (NOM>0) THEN
         IF (.NOT.PROCESS_MESH_NEIGHBORHOOD(NOM)) RETURN  ! If NOM not controlled by current MPI process, abandon search
         OM => MESHES(NOM)
      ELSEIF (IW<=M%N_EXTERNAL_WALL_CELLS .AND. (SF%HT_DIM>1.OR.SF%VARIABLE_THICKNESS)) THEN 
         ! Do not apply HT3D to VARIABLE_THICKNESS exterior boundary
         WRITE(MESSAGE,'(3A)') 'ERROR(437): SURF ',TRIM(SURFACE(WC%SURF_INDEX)%ID),' cannot be applied to an exterior boundary.'
         CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
         RETURN
      ENDIF
   ENDIF

   OLD_THICKNESS = THICKNESS
   SELECT CASE(IOR)
      CASE( 1) ; THICKNESS = ABS(BC%X - OM%X(II))
      CASE(-1) ; THICKNESS = ABS(BC%X - OM%X(II-1))
      CASE( 2) ; THICKNESS = ABS(BC%Y - OM%Y(JJ))
      CASE(-2) ; THICKNESS = ABS(BC%Y - OM%Y(JJ-1))
      CASE( 3) ; THICKNESS = ABS(BC%Z - OM%Z(KK))
      CASE(-3) ; THICKNESS = ABS(BC%Z - OM%Z(KK-1))
   END SELECT

   IC = OM%CELL_INDEX(II,JJ,KK)

   ! For VARIABLE_THICKNESS and HT3D cases, get material information from obstruction

   IF (SF%VARIABLE_THICKNESS .OR. SF%HT_DIM>1) THEN 

      ! Determine the index of the obstruction (OBST_INDEX) that occupies the cell with index IC.
      ! If the obstruction is thin, search all obstructions in the mesh for one that corresponds to the upwind cell face.

      OBST_INDEX_PREVIOUS = OBST_INDEX
      OBST_INDEX = OM%CELL(IC)%OBST_INDEX
      
      IF (ITER==1 .AND. OBST_INDEX<1) THEN
         IIF=II ; JJF=JJ ; KKF=KK
         IF (BC%IOR==-1) IIF=II-1 ; IF (BC%IOR==-2) JJF=JJ-1 ; IF (BC%IOR==-3) KKF=KK-1
         SUCCESS = .FALSE.
         DO OBST_INDEX=1,OM%N_OBST
            OB => OM%OBSTRUCTION(OBST_INDEX)
            IF (OB%I1/=OB%I2 .AND. ABS(BC%IOR)==1) CYCLE
            IF (OB%J1/=OB%J2 .AND. ABS(BC%IOR)==2) CYCLE
            IF (OB%K1/=OB%K2 .AND. ABS(BC%IOR)==3) CYCLE
            IF (IIF>=OB%I1.AND.IIF<=OB%I2.AND.JJF>=OB%J1.AND.JJF<=OB%J2.AND.KKF>=OB%K1.AND.KKF<=OB%K2) THEN
               SUCCESS = .TRUE.
               EXIT
            ENDIF
         ENDDO
         IF (.NOT.SUCCESS) THEN
            WRITE(MESSAGE,'(3A)') 'ERROR(368): SURF ',TRIM(SF%ID),' has a problem with VARIABLE_THICKNESS or HT3D.'
            CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
         ENDIF
      ENDIF

      OB => OM%OBSTRUCTION(OBST_INDEX)

      ! If this is a thin obstruction, use its actual user-specified coordinates to determine THICKNESS

      IF (OB%THIN) THIN_OBSTRUCTION = .TRUE.

      IF (THICKNESS<TWO_EPSILON_EB .AND. OBST_INDEX>0) THEN
         SELECT CASE(ABS(BC%IOR))
            CASE(1) ; IF (OB%I1==OB%I2) THICKNESS = OB%X2 - OB%X1
            CASE(2) ; IF (OB%J1==OB%J2) THICKNESS = OB%Y2 - OB%Y1
            CASE(3) ; IF (OB%K1==OB%K2) THICKNESS = OB%Z2 - OB%Z1
         END SELECT
      ENDIF

      IF (OBST_INDEX>0) THEN
         IF (OB%MATL_INDEX(1)<1) THEN
            IF (ITER==1.AND.SF%N_MATL>0) THEN
               OB%MATL_INDEX(1:SF%N_MATL) = SF%MATL_INDEX(1:SF%N_MATL)
               OB%MATL_MASS_FRACTION(1:SF%N_LAYER_MATL(1)) = SF%MATL_MASS_FRACTION(1,1:SF%N_LAYER_MATL(1))
            ELSE
               WRITE(MESSAGE,'(3A)') 'ERROR(375): OBST ',TRIM(OB%ID),' is VARIABLE_THICKNESS or HT3D and needs a MATL_ID.'
               CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.) ; RETURN
            ENDIF
         ENDIF
         HEAT_SOURCE_OBST(N_LAYERS_OBST) = OB%HEAT_SOURCE
         RAMP_IHS_INDEX_OBST(N_LAYERS_OBST) = OB%RAMP_IHS_INDEX
         IF (OB%STRETCH_FACTOR>0._EB)   STRETCH_FACTOR_OBST(N_LAYERS_OBST)    = OB%STRETCH_FACTOR
         IF (OB%CELL_SIZE>0._EB .AND. .NOT.THIN_OBSTRUCTION) CELL_SIZE_OBST(N_LAYERS_OBST) = OB%CELL_SIZE
         IF (OB%CELL_SIZE_FACTOR>0._EB) CELL_SIZE_FACTOR_OBST(N_LAYERS_OBST)  = OB%CELL_SIZE_FACTOR
         IF (OB%N_LAYER_CELLS_MAX>0)    N_LAYER_CELLS_MAX_OBST(N_LAYERS_OBST) = OB%N_LAYER_CELLS_MAX
      ENDIF

      LAYER_THICKNESS_OBST(N_LAYERS_OBST) = LAYER_THICKNESS_OBST(N_LAYERS_OBST) + THICKNESS - OLD_THICKNESS
      IF (OBST_INDEX>0) CALL ADD_MATERIAL(MAX_MATERIALS,OB%MATL_INDEX,N_MATL_OBST,MATL_INDEX_OBST)
      IF (OBST_INDEX/=OBST_INDEX_PREVIOUS .AND. OBST_INDEX_PREVIOUS>0 .AND. OBST_INDEX>0) THEN
         OB_PREV => OM_PREV%OBSTRUCTION(OBST_INDEX_PREVIOUS)
         IF ( (ANY(OB%MATL_MASS_FRACTION(:)/=OB_PREV%MATL_MASS_FRACTION(:),DIM=1)) .OR. &
              (ANY(OB%MATL_INDEX(:)        /=OB_PREV%MATL_INDEX(:)        ,DIM=1)) ) THEN
            N_LAYERS_OBST = N_LAYERS_OBST + 1
            LAYER_THICKNESS_OBST(N_LAYERS_OBST) = 0._EB
            HEAT_SOURCE_OBST(N_LAYERS_OBST)        = OB%HEAT_SOURCE    
            HEAT_SOURCE_OBST(N_LAYERS_OBST-1)      = OB_PREV%HEAT_SOURCE
            RAMP_IHS_INDEX_OBST(N_LAYERS_OBST) = OB%RAMP_IHS_INDEX
            RAMP_IHS_INDEX_OBST(N_LAYERS_OBST-1) = OB_PREV%RAMP_IHS_INDEX
            IF (OB%STRETCH_FACTOR>0._EB)        STRETCH_FACTOR_OBST(N_LAYERS_OBST)     = OB%STRETCH_FACTOR 
            IF (OB_PREV%STRETCH_FACTOR>0._EB)   STRETCH_FACTOR_OBST(N_LAYERS_OBST-1)   = OB_PREV%STRETCH_FACTOR
            IF (OB%CELL_SIZE>0._EB .AND. .NOT.OB%THIN)      CELL_SIZE_OBST(N_LAYERS_OBST)   = OB%CELL_SIZE
            IF (OB_PREV%CELL_SIZE>0._EB .AND. .NOT.OB%THIN) CELL_SIZE_OBST(N_LAYERS_OBST-1) = OB_PREV%CELL_SIZE
            IF (OB%CELL_SIZE_FACTOR>0._EB)      CELL_SIZE_FACTOR_OBST(N_LAYERS_OBST)   = OB%CELL_SIZE_FACTOR
            IF (OB_PREV%CELL_SIZE_FACTOR>0._EB) CELL_SIZE_FACTOR_OBST(N_LAYERS_OBST-1) = OB_PREV%CELL_SIZE_FACTOR
            IF (OB%N_LAYER_CELLS_MAX>0)         N_LAYER_CELLS_MAX_OBST(N_LAYERS_OBST)  = OB%N_LAYER_CELLS_MAX
            IF (OB_PREV%N_LAYER_CELLS_MAX>0)    N_LAYER_CELLS_MAX_OBST(N_LAYERS_OBST-1)= OB_PREV%N_LAYER_CELLS_MAX
         ENDIF
      ENDIF
 
      IF (OBST_INDEX>0) THEN
         DO NN=1,N_MATL_OBST
            DO NNN=1,MAX_MATERIALS
                  IF (OB%MATL_INDEX(NNN)==MATL_INDEX_OBST(NN)) &
                     MATL_MASS_FRACTION_OBST(N_LAYERS_OBST,NN) = OB%MATL_MASS_FRACTION(NNN)
            ENDDO
         ENDDO
      ENDIF

   ENDIF

   ! Determine if the back face is found

   IF ((.NOT.OM%CELL(IC)%SOLID .AND. OM%CELL(IC)%WALL_INDEX(IOR)>0) .OR. NOM==0) THEN ! the back wall face is found
      IF (NOM>0 .AND. SF%BACKING/=EXPOSED) RETURN  ! No need to assign back cell information for anything but exposed backing
      ONE_D%BACK_INDEX = OM%CELL(IC)%WALL_INDEX(IOR)
      ONE_D%BACK_MESH  = NOM
      ONE_D%BACK_SURF  = OM%CELL(IC)%SURF_INDEX(IOR)
      IF (NOM>0) THEN
         OS => M%OMESH(NOM)%WALL_RECV_BUFFER
         IF (.NOT.ALLOCATED(OS%ITEM_INDEX)) THEN
            OS%N_ITEMS_DIM = 50
            ALLOCATE(OS%ITEM_INDEX(1:OS%N_ITEMS_DIM))
            ALLOCATE(OS%SURF_INDEX(1:OS%N_ITEMS_DIM))
         ENDIF
         IF (COUNT(OS%ITEM_INDEX(1:OS%N_ITEMS)==ONE_D%BACK_INDEX)==0) THEN
            IF (OS%N_ITEMS>=OS%N_ITEMS_DIM) THEN
               CALL REALLOCATE_INTEGER_ARRAY(OS%ITEM_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
               CALL REALLOCATE_INTEGER_ARRAY(OS%SURF_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
               OS%N_ITEMS_DIM = OS%N_ITEMS_DIM + 50
            ENDIF
            OS%N_ITEMS = OS%N_ITEMS + 1
            OS%ITEM_INDEX(OS%N_ITEMS) = ONE_D%BACK_INDEX
            OS%SURF_INDEX(OS%N_ITEMS) = ONE_D%BACK_SURF
         ENDIF
      ENDIF
      EXIT FIND_BACK_WALL_CELL
   ENDIF

   ! If 1-D solid and the user-specified thickness is less than the current thickness, abandon the search for back-wall cell

   IF (.NOT.SF%VARIABLE_THICKNESS .AND. SF%HT_DIM==1 .AND. THICKNESS>SUM(SF%LAYER_THICKNESS)) RETURN

   SELECT CASE(IOR)  ! New cell indices as we march deeper into the obstruction
      CASE(-1) ; II=II+1
      CASE( 1) ; II=II-1
      CASE(-2) ; JJ=JJ+1
      CASE( 2) ; JJ=JJ-1
      CASE(-3) ; KK=KK+1
      CASE( 3) ; KK=KK-1
   END SELECT

ENDDO FIND_BACK_WALL_CELL

! If the user has specified LINING materials (HT3D or VARIABLE_THICKNESS with SURF MATLs and THICKNESS), add this information to 
! existing lists of layers and materials.

IF (SF%VARIABLE_THICKNESS .OR. SF%HT_DIM>1) THEN

   N_LAYERS = 0
   N_MATLS = N_MATL_OBST
   MATL_INDEX(1:N_MATLS) = MATL_INDEX_OBST(1:N_MATLS)  ! MATL_INDEX_OBST is taken from the OBSTs that make up the solid
   MATL_MASS_FRACTION = 0._EB
   LAYER_THICKNESS = 0._EB
   MINIMUM_LAYER_THICKNESS = 0._EB
   HT3D_LAYER = .FALSE.
   FRONT_LINING_THICKNESS = 0._EB
   BACK_LINING_THICKNESS = 0._EB
   IF (SF%N_LAYERS>0 .AND. SF%LINING) THEN
      CALL ADD_MATERIAL(SF%N_MATL,SF%MATL_INDEX,N_MATLS,MATL_INDEX)       ! Add materials from the front surface lining
      IF (SF%LINING) FRONT_LINING_THICKNESS = SUM(SF%LAYER_THICKNESS(1:SF%N_LAYERS))
   ENDIF
   SF_BACK => SURFACE(ONE_D%BACK_SURF)
   IF (SF_BACK%N_LAYERS>0 .AND. SF_BACK%LINING) THEN
      CALL ADD_MATERIAL(SF_BACK%N_MATL,SF_BACK%MATL_INDEX,N_MATLS,MATL_INDEX)  ! Add materials from the back surface lining
      IF (SF_BACK%LINING) BACK_LINING_THICKNESS = SUM(SF_BACK%LAYER_THICKNESS(1:SF_BACK%N_LAYERS))
   ENDIF

   ! Offset the coordinates of the wall face to account for linings

   IF (THIN_OBSTRUCTION) THEN
      SELECT CASE(BC%IOR)
         CASE(-1) ; BC%X1 = BC%X1 - FRONT_LINING_THICKNESS
         CASE( 1) ; BC%X1 = BC%X2 + FRONT_LINING_THICKNESS
         CASE(-2) ; BC%Y1 = BC%Y1 - FRONT_LINING_THICKNESS
         CASE( 2) ; BC%Y1 = BC%Y2 + FRONT_LINING_THICKNESS
         CASE(-3) ; BC%Z1 = BC%Z1 - FRONT_LINING_THICKNESS
         CASE( 3) ; BC%Z1 = BC%Z2 + FRONT_LINING_THICKNESS
      END SELECT
   ENDIF

   ! Copy the front face SURF layer information into the holding arrays

   DO NL=1,SF%N_LAYERS
      IF (.NOT.SF%LINING) EXIT
      N_LAYERS = N_LAYERS + 1
      LAYER_THICKNESS(N_LAYERS)  = SF%LAYER_THICKNESS(N_LAYERS) 
      MINIMUM_LAYER_THICKNESS(N_LAYERS)  = SF%MINIMUM_LAYER_THICKNESS(N_LAYERS) 
      HT3D_LAYER(N_LAYERS)       = .FALSE.
      HEAT_SOURCE(N_LAYERS)      = SF%HEAT_SOURCE(N_LAYERS)
      RAMP_IHS_INDEX(N_LAYERS)   = SF%RAMP_IHS_INDEX(N_LAYERS)
      STRETCH_FACTOR(N_LAYERS)   = SF%STRETCH_FACTOR(N_LAYERS)
      CELL_SIZE(N_LAYERS)        = SF%CELL_SIZE(N_LAYERS)
      CELL_SIZE_FACTOR(N_LAYERS) = SF%CELL_SIZE_FACTOR(N_LAYERS)
      N_LAYER_CELLS_MAX(N_LAYERS) = SF%N_LAYER_CELLS_MAX(N_LAYERS)
      DO NN=1,SF%N_LAYER_MATL(NL)
         DO NNN=1,N_MATLS
            IF (SF%LAYER_MATL_INDEX(NL,NN)==MATL_INDEX(NNN)) MATL_MASS_FRACTION(N_LAYERS,NNN) = SF%MATL_MASS_FRACTION(NL,NN)
         ENDDO
      ENDDO
   ENDDO

   ! Add layers that are associated with the underlying OBSTructions

   IF (.NOT.THIN_OBSTRUCTION) THEN
      LAYER_THICKNESS_OBST_TOTAL = SUM(LAYER_THICKNESS_OBST(1:N_LAYERS_OBST))
      LAYER_THICKNESS_OBST(1:N_LAYERS_OBST) = LAYER_THICKNESS_OBST(1:N_LAYERS_OBST)*&
         (LAYER_THICKNESS_OBST_TOTAL-FRONT_LINING_THICKNESS-BACK_LINING_THICKNESS)/LAYER_THICKNESS_OBST_TOTAL
   ENDIF

   DO NL=1,N_LAYERS_OBST
      N_LAYERS = N_LAYERS + 1
      LAYER_THICKNESS(N_LAYERS)  = LAYER_THICKNESS_OBST(NL)
      MINIMUM_LAYER_THICKNESS(N_LAYERS) = SF%MINIMUM_LAYER_THICKNESS(1)
      HT3D_LAYER(N_LAYERS)       = .TRUE.
      HEAT_SOURCE(N_LAYERS)      = HEAT_SOURCE_OBST(NL)
      RAMP_IHS_INDEX(N_LAYERS)   = RAMP_IHS_INDEX_OBST(NL)
      STRETCH_FACTOR(N_LAYERS)   = STRETCH_FACTOR_OBST(NL)
      CELL_SIZE(N_LAYERS)        = CELL_SIZE_OBST(NL)
      CELL_SIZE_FACTOR(N_LAYERS) = CELL_SIZE_FACTOR_OBST(NL)
      N_LAYER_CELLS_MAX(N_LAYERS) = N_LAYER_CELLS_MAX_OBST(NL)
      DO NN=1,N_MATL_OBST
         DO NNN=1,N_MATLS
            IF (MATL_INDEX_OBST(NN)==MATL_INDEX(NNN)) MATL_MASS_FRACTION(N_LAYERS,NNN) = MATL_MASS_FRACTION_OBST(NL,NN)
         ENDDO
      ENDDO
   ENDDO

   ! Add layers from the back surface lining

   DO NL=1,SF_BACK%N_LAYERS
      IF (.NOT.SF_BACK%LINING) EXIT
      N_LAYERS = N_LAYERS + 1
      LAYER_THICKNESS(N_LAYERS)  = SF_BACK%LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1)
      MINIMUM_LAYER_THICKNESS(N_LAYERS)  = SF_BACK%MINIMUM_LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1)
      HT3D_LAYER(N_LAYERS)       = .FALSE.
      HEAT_SOURCE(N_LAYERS)      = SF_BACK%HEAT_SOURCE(SF_BACK%N_LAYERS-NL+1)
      RAMP_IHS_INDEX(N_LAYERS)   = SF_BACK%RAMP_IHS_INDEX(SF_BACK%N_LAYERS-NL+1)
      STRETCH_FACTOR(N_LAYERS)   = SF_BACK%STRETCH_FACTOR(SF_BACK%N_LAYERS-NL+1)
      CELL_SIZE(N_LAYERS)        = SF_BACK%CELL_SIZE(SF_BACK%N_LAYERS-NL+1)
      CELL_SIZE_FACTOR(N_LAYERS) = SF_BACK%CELL_SIZE_FACTOR(SF_BACK%N_LAYERS-NL+1)
      N_LAYER_CELLS_MAX(N_LAYERS)= SF_BACK%N_LAYER_CELLS_MAX(SF_BACK%N_LAYERS-NL+1)
      DO NN=1,SF_BACK%N_LAYER_MATL(NL)
         DO NNN=1,N_MATLS
            IF (SF_BACK%LAYER_MATL_INDEX(SF_BACK%N_LAYERS-NL+1,NN)==MATL_INDEX(NNN)) &
               MATL_MASS_FRACTION(N_LAYERS,NNN) = SF_BACK%MATL_MASS_FRACTION(SF_BACK%N_LAYERS-NL+1,NN)
         ENDDO
      ENDDO
   ENDDO

   ! Reallocate ONE_D arrays holding layer and material info for HT3D and VARIABLE_THICKNESS objects

   ONE_D%N_LAYERS = N_LAYERS
   ONE_D%N_MATL   = N_MATLS
   DEALLOCATE(ONE_D%MATL_COMP)       ; ALLOCATE(ONE_D%MATL_COMP(ONE_D%N_MATL))
   DEALLOCATE(ONE_D%MATL_INDEX)      ; ALLOCATE(ONE_D%MATL_INDEX(ONE_D%N_MATL))
   DEALLOCATE(ONE_D%LAYER_THICKNESS) ; ALLOCATE(ONE_D%LAYER_THICKNESS(ONE_D%N_LAYERS))
   DEALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS) ; ALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS(ONE_D%N_LAYERS))
   DEALLOCATE(ONE_D%HT3D_LAYER)      ; ALLOCATE(ONE_D%HT3D_LAYER(ONE_D%N_LAYERS))
   ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS) = LAYER_THICKNESS(1:ONE_D%N_LAYERS)
   ONE_D%MINIMUM_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = MINIMUM_LAYER_THICKNESS(1:ONE_D%N_LAYERS)
   ONE_D%HT3D_LAYER(1:ONE_D%N_LAYERS) = HT3D_LAYER(1:ONE_D%N_LAYERS)
   DO NN=1,ONE_D%N_MATL
      ALLOCATE(ONE_D%MATL_COMP(NN)%MASS_FRACTION(ONE_D%N_LAYERS))
      ONE_D%MATL_INDEX(NN) = MATL_INDEX(NN)
      IF (MATERIAL(ONE_D%MATL_INDEX(NN))%PYROLYSIS_MODEL/=PYROLYSIS_NONE) THEN
         ONE_D%PYROLYSIS_MODEL = PYROLYSIS_PREDICTED
         SF%SPECIES_BC_INDEX = SPECIFIED_MASS_FLUX
      ENDIF
      DO NL=1,ONE_D%N_LAYERS
         ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL) = MATL_MASS_FRACTION(NL,NN)
      ENDDO
      IF (MATERIAL(ONE_D%MATL_INDEX(NN))%KAPPA_S<4.9E4_EB) ONE_D%INTERNAL_RADIATION = .TRUE.
   ENDDO
   DEALLOCATE(ONE_D%HEAT_SOURCE)       ; ALLOCATE(ONE_D%HEAT_SOURCE(ONE_D%N_LAYERS))       ; ONE_D%HEAT_SOURCE = 0._EB
   DEALLOCATE(ONE_D%RAMP_IHS_INDEX)    ; ALLOCATE(ONE_D%RAMP_IHS_INDEX(ONE_D%N_LAYERS))    ; ONE_D%RAMP_IHS_INDEX = 0._EB
   DEALLOCATE(ONE_D%STRETCH_FACTOR)    ; ALLOCATE(ONE_D%STRETCH_FACTOR(ONE_D%N_LAYERS))    ; ONE_D%STRETCH_FACTOR =2._EB
   DEALLOCATE(ONE_D%CELL_SIZE)         ; ALLOCATE(ONE_D%CELL_SIZE(ONE_D%N_LAYERS))         ; ONE_D%CELL_SIZE = -1._EB
   DEALLOCATE(ONE_D%CELL_SIZE_FACTOR)  ; ALLOCATE(ONE_D%CELL_SIZE_FACTOR(ONE_D%N_LAYERS))  ; ONE_D%CELL_SIZE_FACTOR = 1._EB
   DEALLOCATE(ONE_D%N_LAYER_CELLS_MAX) ; ALLOCATE(ONE_D%N_LAYER_CELLS_MAX(ONE_D%N_LAYERS)) ; ONE_D%N_LAYER_CELLS_MAX = 999
   ONE_D%HEAT_SOURCE(1:ONE_D%N_LAYERS)       = HEAT_SOURCE(1:ONE_D%N_LAYERS)
   ONE_D%RAMP_IHS_INDEX(1:ONE_D%N_LAYERS)    = RAMP_IHS_INDEX(1:ONE_D%N_LAYERS)
   ONE_D%STRETCH_FACTOR(1:ONE_D%N_LAYERS)    = STRETCH_FACTOR(1:ONE_D%N_LAYERS)
   ONE_D%CELL_SIZE(1:ONE_D%N_LAYERS)         = CELL_SIZE(1:ONE_D%N_LAYERS)
   ONE_D%CELL_SIZE_FACTOR(1:ONE_D%N_LAYERS)  = CELL_SIZE_FACTOR(1:ONE_D%N_LAYERS)
   ONE_D%N_LAYER_CELLS_MAX(1:ONE_D%N_LAYERS) = N_LAYER_CELLS_MAX(1:ONE_D%N_LAYERS)

ENDIF

END SUBROUTINE FIND_WALL_BACK_INDEX


!> \brief Find back index of thin wall
!> \param NM Mesh number
!> \param ITW Thin wall index
!> \details ITW is the index of a thin wall cell, which can be thought of as segment of the edge of a single thin obstruction. 
!> This routine marches from edge to opposite edge looking for the "back" thin wall index.

SUBROUTINE FIND_THIN_WALL_BACK_INDEX(NM,ITW)

USE GEOMETRY_FUNCTIONS, ONLY: SEARCH_OTHER_MESHES
USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_INTEGER_ARRAY
INTEGER, INTENT(IN) :: NM,ITW
INTEGER :: II,JJ,KK,IC,IOR,IEC,ITW2,NOM,NN,NNN,NL,N_MATLS,IIGM,IIGP,JJGM,JJGP,KKGM,KKGP,ICM,ICP
INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX
REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION
REAL(EB) :: XXC,YYC,ZZC
TYPE (MESH_TYPE), POINTER :: M
TYPE (THIN_WALL_TYPE), POINTER :: TW
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE (MESH_TYPE), POINTER :: OM
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (STORAGE_TYPE), POINTER :: OS

M => MESHES(NM)

TW => M%THIN_WALL(ITW)
BC => M%BOUNDARY_COORD(TW%BC_INDEX)
OB => M%OBSTRUCTION(TW%OBST_INDEX)
SF => SURFACE(TW%SURF_INDEX)
N_MATLS = 0
MATL_INDEX = 0

! If OBSTstruction to which the thin wall cell is attached has no material index, use the indices associated with the SURF

IF (OB%MATL_INDEX(1)<1) THEN
   OB%MATL_INDEX(1:SF%N_MATL) = SF%MATL_INDEX(1:SF%N_MATL)
   OB%MATL_MASS_FRACTION(1:SF%N_LAYER_MATL(1)) = SF%MATL_MASS_FRACTION(1,1:SF%N_LAYER_MATL(1))
ENDIF

! Form an array of N_MATLS material indices, MATL_INDEX, for this thin wall cell. This
! list accounts for all materials associated with the OBSTs and SURFs along the distance through the solid.

CALL ADD_MATERIAL(MAX_MATERIALS,OB%MATL_INDEX,N_MATLS,MATL_INDEX)

! A thin wall cell only has one layer and one obstruction. This loop transfers the material mass fractions 
! from the OBST to the save array.

DO NN=1,N_MATLS
   DO NNN=1,MAX_MATERIALS
         IF (OB%MATL_INDEX(NNN)==MATL_INDEX(NN)) MATL_MASS_FRACTION(1,NN) = OB%MATL_MASS_FRACTION(NNN)
   ENDDO
ENDDO

II  = BC%II
JJ  = BC%JJ
KK  = BC%KK
IOR = BC%IOR
IEC = TW%IEC
NOM = NM
OM => MESHES(NOM)

! Find one or two WALL cells that are adjacent to this THIN_WALL cell

IIGM=II ; JJGM=JJ ; KKGM = KK ; IIGP=II ; JJGP=JJ ; KKGP = KK
SELECT CASE(IEC)
   CASE(1)
      SELECT CASE(IOR)
         CASE(-2) ; JJGM=JJ   ; JJGP=JJ   ; KKGM=KK   ; KKGP=KK+1
         CASE( 2) ; JJGM=JJ+1 ; JJGP=JJ+1 ; KKGM=KK   ; KKGP=KK+1
         CASE(-3) ; JJGM=JJ   ; JJGP=JJ+1 ; KKGM=KK   ; KKGP=KK
         CASE( 3) ; JJGM=JJ   ; JJGP=JJ+1 ; KKGM=KK+1 ; KKGP=KK+1
      END SELECT
   CASE(2)
      SELECT CASE(IOR)
         CASE(-1) ; IIGM=II   ; IIGP=II   ; KKGM=KK   ; KKGP=KK+1
         CASE( 1) ; IIGM=II+1 ; IIGP=II+1 ; KKGM=KK   ; KKGP=KK+1
         CASE(-3) ; IIGM=II   ; IIGP=II+1 ; KKGM=KK   ; KKGP=KK
         CASE( 3) ; IIGM=II   ; IIGP=II+1 ; KKGM=KK+1 ; KKGP=KK+1
      END SELECT
   CASE(3)
      SELECT CASE(IOR)
         CASE(-1) ; IIGM=II   ; IIGP=II   ; JJGM=JJ   ; JJGP=JJ+1
         CASE( 1) ; IIGM=II+1 ; IIGP=II+1 ; JJGM=JJ   ; JJGP=JJ+1
         CASE(-2) ; IIGM=II   ; IIGP=II+1 ; JJGM=JJ   ; JJGP=JJ
         CASE( 2) ; IIGM=II   ; IIGP=II+1 ; JJGM=JJ+1 ; JJGP=JJ+1
      END SELECT
END SELECT

ICM = M%CELL_INDEX(IIGM,JJGM,KKGM)
ICP = M%CELL_INDEX(IIGP,JJGP,KKGP)
TW%WALL_INDEX_M = M%CELL(ICM)%WALL_INDEX(-IOR)
TW%WALL_INDEX_P = M%CELL(ICP)%WALL_INDEX(-IOR)

! Look for the back THIN_WALL cell; that is, the thin wall cell on the other side of the obstruction

FIND_BACK_THIN_WALL_CELL: DO

   IF ((II==0.AND.IOR==1) .OR. (II==OM%IBAR.AND.IOR==-1) .OR. &
       (JJ==0.AND.IOR==2) .OR. (JJ==OM%JBAR.AND.IOR==-2) .OR. &
       (KK==0.AND.IOR==3) .OR. (KK==OM%KBAR.AND.IOR==-3)) THEN
      XXC=OM%XC(II) ; YYC=OM%YC(JJ) ; ZZC=OM%ZC(KK)
      IF (II==0      .AND.IOR== 1) XXC = OM%X(II) - MESH_SEPARATION_DISTANCE
      IF (II==OM%IBAR.AND.IOR==-1) XXC = OM%X(II) + MESH_SEPARATION_DISTANCE
      IF (JJ==0      .AND.IOR== 2) YYC = OM%Y(JJ) - MESH_SEPARATION_DISTANCE
      IF (JJ==OM%JBAR.AND.IOR==-2) YYC = OM%Y(JJ) + MESH_SEPARATION_DISTANCE
      IF (KK==0      .AND.IOR== 3) ZZC = OM%Z(KK) - MESH_SEPARATION_DISTANCE
      IF (KK==OM%KBAR.AND.IOR==-3) ZZC = OM%Z(KK) + MESH_SEPARATION_DISTANCE
      CALL SEARCH_OTHER_MESHES(XXC,YYC,ZZC,NOM,II,JJ,KK)
      IF (NOM==0) RETURN
      OM => MESHES(NOM)
   ENDIF

   ! Look for the other side of the thin obstruction, using its ORDINAL value as a unique identifier

   IC = OM%CELL_INDEX(II,JJ,KK)
   ITW2 = OM%CELL(IC)%THIN_WALL_INDEX(-IOR,IEC)

   IF (ITW2>0 .AND. OM%CELL(IC)%THIN_OBST_INDEX(-IOR,IEC)==OB%ORDINAL) THEN ! the back thin wall is found
      ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX)
      ONE_D%BACK_INDEX = OM%CELL(IC)%THIN_WALL_INDEX(-IOR,IEC)
      ONE_D%BACK_MESH  = NOM
      ONE_D%BACK_SURF  = TW%SURF_INDEX
      OS => M%OMESH(NOM)%THIN_WALL_RECV_BUFFER
      IF (.NOT.ALLOCATED(OS%ITEM_INDEX)) THEN
         OS%N_ITEMS_DIM = 50
         ALLOCATE(OS%ITEM_INDEX(1:OS%N_ITEMS_DIM))
         ALLOCATE(OS%SURF_INDEX(1:OS%N_ITEMS_DIM))
      ENDIF
      IF (COUNT(OS%ITEM_INDEX(1:OS%N_ITEMS)==ONE_D%BACK_INDEX)==0) THEN
         IF (OS%N_ITEMS>=OS%N_ITEMS_DIM) THEN
            CALL REALLOCATE_INTEGER_ARRAY(OS%ITEM_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
            CALL REALLOCATE_INTEGER_ARRAY(OS%SURF_INDEX,1,OS%N_ITEMS_DIM,OS%N_ITEMS_DIM+50)
            OS%N_ITEMS_DIM = OS%N_ITEMS_DIM + 50
         ENDIF
         OS%N_ITEMS = OS%N_ITEMS + 1
         OS%ITEM_INDEX(OS%N_ITEMS) = ONE_D%BACK_INDEX
         OS%SURF_INDEX(OS%N_ITEMS) = ONE_D%BACK_SURF
      ENDIF
      SELECT CASE(ABS(IOR))
         CASE(1) ; ONE_D%LAYER_THICKNESS(1) = OB%UNDIVIDED_INPUT_LENGTH(1)
         CASE(2) ; ONE_D%LAYER_THICKNESS(1) = OB%UNDIVIDED_INPUT_LENGTH(2)
         CASE(3) ; ONE_D%LAYER_THICKNESS(1) = OB%UNDIVIDED_INPUT_LENGTH(3)
      END SELECT
      IF (OB%CELL_SIZE>0._EB) THEN
         ONE_D%CELL_SIZE(1) = OB%CELL_SIZE
         ONE_D%STRETCH_FACTOR(1) = 1._EB
      ENDIF
      EXIT FIND_BACK_THIN_WALL_CELL
   ENDIF

   ! If the back thin wall index is not found, update the cell indices and continue marching deeper into the obstruction

   SELECT CASE(IOR)
      CASE(-1) ; II=II+1
      CASE( 1) ; II=II-1
      CASE(-2) ; JJ=JJ+1
      CASE( 2) ; JJ=JJ-1
      CASE(-3) ; KK=KK+1
      CASE( 3) ; KK=KK-1
   END SELECT

ENDDO FIND_BACK_THIN_WALL_CELL

! Take the array of MATL_INDEX and MATL_MASS_FRACTION and save them in the ONE_D derived type variable.

ONE_D%N_MATL = N_MATLS
DEALLOCATE(ONE_D%MATL_COMP) ; ALLOCATE(ONE_D%MATL_COMP(ONE_D%N_MATL))
DEALLOCATE(ONE_D%MATL_INDEX) ; ALLOCATE(ONE_D%MATL_INDEX(ONE_D%N_MATL))
DO NN=1,ONE_D%N_MATL
   ALLOCATE(ONE_D%MATL_COMP(NN)%MASS_FRACTION(ONE_D%N_LAYERS))
   ONE_D%MATL_INDEX(NN) = MATL_INDEX(NN)
   DO NL=1,ONE_D%N_LAYERS
      ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL) = MATL_MASS_FRACTION(NL,NN)
   ENDDO
ENDDO

END SUBROUTINE FIND_THIN_WALL_BACK_INDEX


!> \brief Update list of material indices
!> \details The list of materials on the search list are checked against the X list and added if not there. Then the residues
!> of the materials added to the X list are checked, and the residues of the residues, etc.
!> \param N_MATLS_SEARCH Number of materials in the array to be searched
!> \param MATL_INDEX_SEARCH Array of material indices on the search list
!> \param N_MATLS_X Number of materials on the new list
!> \param MATL_INDEX_X Array of new material indices

SUBROUTINE ADD_MATERIAL(N_MATLS_SEARCH,MATL_INDEX_SEARCH,N_MATLS_X,MATL_INDEX_X)

INTEGER, INTENT(IN) :: N_MATLS_SEARCH
INTEGER, INTENT(IN), DIMENSION(N_MATLS_SEARCH) :: MATL_INDEX_SEARCH
INTEGER, INTENT(INOUT) :: N_MATLS_X
INTEGER, INTENT(INOUT), DIMENSION(MAX_MATERIALS) :: MATL_INDEX_X
INTEGER :: JJ,MATL_SUM(N_MATL)

MATL_SUM = 0
MATL_LOOP: DO JJ=1,N_MATLS_SEARCH
    IF (MATL_INDEX_SEARCH(JJ) < 1) EXIT MATL_LOOP
    MATL_SUM = MATL_SUM + MATERIAL(MATL_INDEX_SEARCH(JJ))%CHILD_MATL
ENDDO MATL_LOOP

WHERE (MATL_SUM > 1) MATL_SUM = 1

MATL_LOOP_2: DO JJ=1,N_MATL
   IF (MATL_SUM(JJ) == 1) THEN
      IF (FINDLOC(MATL_INDEX_X,JJ,DIM=1)==0) THEN
         N_MATLS_X = N_MATLS_X + 1
         MATL_INDEX_X(N_MATLS_X) = JJ
      ENDIF
   ENDIF
ENDDO MATL_LOOP_2

RETURN

END SUBROUTINE ADD_MATERIAL


!> \brief Check to see if a cell or OBSTruction is to be created or removed, or a VENT activated of deactivated
!> \param T Current time (s)
!> \param NM Mesh number

SUBROUTINE OPEN_AND_CLOSE(T,NM)

USE MESH_POINTERS
USE MEMORY_FUNCTIONS, ONLY : RE_ALLOCATE_STRINGS
USE CONTROL_VARIABLES, ONLY : CONTROL
USE DEVICE_VARIABLES, ONLY : DEVICE
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
REAL(EB), INTENT(IN) :: T
REAL(EB) :: TNOW
INTEGER  :: N,II,JJ,KK,IW,IC,VENT_INDEX,CVENT_INDEX
INTEGER, INTENT(IN) :: NM
LOGICAL :: CREATE_OBST,REMOVE_OBST,ACTIVATE_VENT,DEACTIVATE_VENT,ANY_REMOVE_OBST
CHARACTER(12) :: SV_LABEL
TYPE (VENTS_TYPE), POINTER :: VT
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1

TNOW = CURRENT_TIME()

CALL POINT_TO_MESH(NM)

ANY_REMOVE_OBST = .FALSE.

! Check to see if an obstacle is to be removed or created

OBST_LOOP: DO N=1,N_OBST

   OB=>OBSTRUCTION(N)
   IF (.NOT. OB%REMOVABLE) CYCLE OBST_LOOP
   CREATE_OBST = .FALSE.
   REMOVE_OBST = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   CREATE_REMOVE_IF:IF (OB%CONSUMABLE .AND. OB%MASS<TWO_EPSILON_EB) THEN
      REMOVE_OBST = .TRUE.
   ELSE CREATE_REMOVE_IF
      SET_T_BEGIN_IF: IF (T<=T_BEGIN) THEN
         ! Set initial state of OBST
         HOLE_FILL_IF: IF (.NOT. OB%HOLE_FILLER) THEN
            !OBST is not a HOLE
            IF (OB%DEVC_INDEX > 0) THEN
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ELSEIF (OB%CTRL_INDEX > 0) THEN
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ENDIF
         ELSE HOLE_FILL_IF
            !OBST is a HOLE.  CREATE/REMOVE also depends on parent OBST.
            CHECK_PARENT: IF (OB%DEVC_INDEX_O > 0 .OR. OB%CTRL_INDEX_O > 0) THEN
               !Parent OBST controllable, check state and if parent OBST is hidden, do not fill hole.
               IF (OB%DEVC_INDEX_O > 0) THEN
                  IF (.NOT. DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ELSEIF(OB%CTRL_INDEX_O > 0) THEN
                  IF (.NOT. CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ENDIF
               !If parent OBST is visible, check to see if hole needs to be made.
               IF (.NOT. REMOVE_OBST) THEN
                  IF (OB%DEVC_INDEX > 0) THEN
                     IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                        REMOVE_OBST = .TRUE.
                     ELSE
                        CREATE_OBST = .TRUE.
                     ENDIF
                  ELSEIF (OB%CTRL_INDEX > 0) THEN
                     IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                        REMOVE_OBST = .TRUE.
                     ELSE
                        CREATE_OBST = .TRUE.
                     ENDIF
                  ENDIF
               ENDIF
            ELSE CHECK_PARENT
               !Parent OBST always present
               IF (OB%DEVC_INDEX > 0) THEN
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ELSEIF (OB%CTRL_INDEX > 0) THEN
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ENDIF
            ENDIF CHECK_PARENT

         ENDIF HOLE_FILL_IF
      ELSE SET_T_BEGIN_IF
      ! Decide if a DEVICE/CONTROL action is needed
         HOLE_FILL_IF_2: IF (.NOT. OB%HOLE_FILLER) THEN
            !OBST is not a HOLE
            IF (OB%DEVC_INDEX > 0) THEN
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
               IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ELSEIF (OB%CTRL_INDEX > 0) THEN
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
               IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                  CREATE_OBST = .TRUE.
               ELSE
                  REMOVE_OBST = .TRUE.
               ENDIF
            ENDIF
         ELSE HOLE_FILL_IF_2
            !OBST is a HOLE.  CREATE/REMOVE also depends on parent OBST.
            CHECK_PARENT_2: IF (OB%DEVC_INDEX_O > 0 .OR. OB%CTRL_INDEX_O > 0) THEN
               !Parent OBST controllable, check state and if parent OBST is hidden, do not fill hole.
               IF (OB%DEVC_INDEX_O > 0) THEN
                  IF (.NOT. DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ELSEIF(OB%CTRL_INDEX_O > 0) THEN
                  IF (.NOT. CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE) REMOVE_OBST = .TRUE.
               ENDIF
               !If parent OBST is visible, check to see if hole needs to be made.
               IF (.NOT. REMOVE_OBST) THEN
                  IF (OB%DEVC_INDEX > 0) THEN
                     IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) THEN
                        IF (OB%DEVC_INDEX_O > 0 .AND. .NOT. DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           IF (DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE .NEQV. DEVICE(OB%DEVC_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ELSEIF(OB%CTRL_INDEX_O > 0 .AND. .NOT. DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           IF (CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE .NEQV. CONTROL(OB%CTRL_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ENDIF
                     ELSE
                        IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                           REMOVE_OBST = .TRUE.
                        ELSE
                           CREATE_OBST = .TRUE.
                        ENDIF
                     ENDIF
                  ELSEIF (OB%CTRL_INDEX > 0) THEN
                     IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) THEN
                        IF (OB%DEVC_INDEX_O > 0 .AND. .NOT. CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           IF (DEVICE(OB%DEVC_INDEX_O)%CURRENT_STATE .NEQV. DEVICE(OB%DEVC_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ELSEIF(OB%CTRL_INDEX_O > 0 .AND. .NOT. CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           IF (CONTROL(OB%CTRL_INDEX_O)%CURRENT_STATE .NEQV. CONTROL(OB%CTRL_INDEX_O)%PRIOR_STATE) &
                               CREATE_OBST=.TRUE.
                        ENDIF
                     ELSE
                        IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                           REMOVE_OBST = .TRUE.
                        ELSE
                           CREATE_OBST = .TRUE.
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ELSE CHECK_PARENT_2
               !Parent OBST not controllable and is always present
               IF (OB%DEVC_INDEX > 0) THEN
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(OB%DEVC_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
                  IF (DEVICE(OB%DEVC_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ELSEIF (OB%CTRL_INDEX > 0) THEN
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(OB%CTRL_INDEX)%PRIOR_STATE) CYCLE OBST_LOOP
                  IF (CONTROL(OB%CTRL_INDEX)%CURRENT_STATE) THEN
                     REMOVE_OBST = .TRUE.
                  ELSE
                     CREATE_OBST = .TRUE.
                  ENDIF
               ENDIF
            ENDIF CHECK_PARENT_2

         ENDIF HOLE_FILL_IF_2

      ENDIF SET_T_BEGIN_IF
   ENDIF CREATE_REMOVE_IF

   SV_LABEL  = 'null'

   IF (CREATE_OBST .AND. OB%HIDDEN) THEN
      OB%HIDDEN = .FALSE.
      SV_LABEL  = 'SHOW_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,1,N)
   ENDIF

   IF (REMOVE_OBST .AND. (.NOT. OB%HIDDEN)) THEN
      OB%HIDDEN = .TRUE.
      SV_LABEL  = 'HIDE_OBST'
      CALL CREATE_OR_REMOVE_OBST(NM,OB%I1,OB%I2,OB%J1,OB%J2,OB%K1,OB%K2,0,N)
   ENDIF

   ! Write a message to the Smokeview .smv file that the obstruction has been created or removed

   IF (SV_LABEL /= 'null') THEN
      IF (N_STRINGS+2>N_STRINGS_MAX) THEN
         CALL RE_ALLOCATE_STRINGS(NM)
         STRING => MESHES(NM)%STRING
      ENDIF
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
      N_STRINGS = N_STRINGS + 1
      WRITE(STRING(N_STRINGS),'(I6,F14.6)') N,T_BEGIN+(T-T_BEGIN)*TIME_SHRINK_FACTOR
   ENDIF

   ! If any REMOVE_OBST store

   IF (REMOVE_OBST) ANY_REMOVE_OBST=.TRUE.

ENDDO OBST_LOOP

! Check to see if a vent should be activated or deactivated

VENT_INDEX  = 0
CVENT_INDEX = 0

VENT_LOOP: DO N=1,N_VENT
   VT => VENTS(N)

   IF (VT%RADIUS>0._EB) THEN
      CVENT_INDEX = CVENT_INDEX + 1
   ELSE
      VENT_INDEX = VENT_INDEX + 1
   ENDIF

   ACTIVATE_VENT   = .FALSE.
   DEACTIVATE_VENT = .FALSE.

   ! Over-ride DEVICE/CONTROL logic

   IF (.NOT.VT%ACTIVATED .AND. T<=T_BEGIN) DEACTIVATE_VENT = .TRUE.

   ! If the VENT is tied to a specific OBST, and the OBST is HIDDEN (not HIDDEN), and the VENT is activated (not activated),
   ! deactivate (activate) the vent.

   IF (VT%OBST_INDEX>0 .AND. OBSTRUCTION(VT%OBST_INDEX)%HIDDEN .AND. VT%ACTIVATED) THEN
      VT%ACTIVATED = .FALSE.
      DEACTIVATE_VENT = .TRUE.
   ENDIF

   IF (VT%OBST_INDEX>0 .AND. .NOT.OBSTRUCTION(VT%OBST_INDEX)%HIDDEN .AND. .NOT.VT%ACTIVATED) THEN
      VT%ACTIVATED = .TRUE.
      ACTIVATE_VENT = .TRUE.
   ENDIF

   ! Decide if a VENT is to activate or de-activate based on a DEVICE or CONTROLLER

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) THEN
      IF (VT%DEVC_INDEX > 0) THEN
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE .EQV. DEVICE(VT%DEVC_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (DEVICE(VT%DEVC_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ELSEIF (VT%CTRL_INDEX > 0) THEN
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE .EQV. CONTROL(VT%CTRL_INDEX)%PRIOR_STATE) CYCLE VENT_LOOP
         IF (CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) THEN
            ACTIVATE_VENT   = .TRUE.
         ELSE
            DEACTIVATE_VENT = .TRUE.
         ENDIF
      ENDIF
   ENDIF

   IF (.NOT.ACTIVATE_VENT .AND. .NOT.DEACTIVATE_VENT) CYCLE VENT_LOOP

   ! Find the wall indices (IW) for the vent and set the activation time (B1%T_IGN) for each one

   DO KK=VT%K1+1,MAX(VT%K1+1,VT%K2)
      DO JJ=VT%J1+1,MAX(VT%J1+1,VT%J2)
         DO II=VT%I1+1,MAX(VT%I1+1,VT%I2)
            SELECT CASE(VT%IOR)
               CASE(1:)
                  IC = CELL_INDEX(II,JJ,KK)
               CASE(-1)
                  IC = CELL_INDEX(II-1,JJ,KK)
               CASE(-2)
                  IC = CELL_INDEX(II,JJ-1,KK)
               CASE(-3)
                  IC = CELL_INDEX(II,JJ,KK-1)
            END SELECT
            IW = CELL(IC)%WALL_INDEX(-VT%IOR)
            IF (IW==0) CYCLE
            B1 => MESHES(NM)%BOUNDARY_PROP1(MESHES(NM)%WALL(IW)%B1_INDEX)
            BC => MESHES(NM)%BOUNDARY_COORD(MESHES(NM)%WALL(IW)%BC_INDEX)

            IF (ACTIVATE_VENT) THEN
               IF (VT%FIRE_SPREAD_RATE>0._EB) THEN
                  B1%T_IGN = T + SQRT((BC%X-VT%X0)**2 + (BC%Y-VT%Y0)**2 + (BC%Z-VT%Z0)**2)/VT%FIRE_SPREAD_RATE
               ELSE
                  B1%T_IGN = T
               ENDIF
            ELSE
               B1%T_IGN = 1.E6_EB
            ENDIF
         ENDDO
      ENDDO
   ENDDO

   ! Write message to .smv file

   IF (VT%RADIUS<0._EB) THEN
      IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_VENT'
      IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_VENT'
   ELSE
      IF (ACTIVATE_VENT)   SV_LABEL = 'OPEN_CVENT'
      IF (DEACTIVATE_VENT) SV_LABEL = 'CLOSE_CVENT'
   ENDIF

   IF (N_STRINGS+2>N_STRINGS_MAX) THEN
      CALL RE_ALLOCATE_STRINGS(NM)
      STRING => MESHES(NM)%STRING
   ENDIF
   N_STRINGS = N_STRINGS + 1
   WRITE(STRING(N_STRINGS),'(A,I3)') SV_LABEL,NM
   N_STRINGS = N_STRINGS + 1
   IF (VT%RADIUS>0._EB) WRITE(STRING(N_STRINGS),'(I6,F10.2)') CVENT_INDEX,T
   IF (VT%RADIUS<0._EB) WRITE(STRING(N_STRINGS),'(I6,F10.2)')  VENT_INDEX,T

ENDDO VENT_LOOP

T_USED(6) = T_USED(6) + CURRENT_TIME() - TNOW
END SUBROUTINE OPEN_AND_CLOSE


!> Create or remove the obstruction whose NODES (not cells) are given by I1, I2, etc.
!> \param NM Mesh number
!> \param I1 Lower x-index of obstruction
!> \param I2 Upper x-index of obstruction
!> \param J1 Lower y-index of obstruction
!> \param J2 Upper y-index of obstruction
!> \param K1 Lower z-index of obstruction
!> \param K2 Upper z-index of obstruction
!> \param CR_INDEX 1 if obstruction is to be created; 0 if removed
!> \param OBST_INDEX Index of the obstruction

SUBROUTINE CREATE_OR_REMOVE_OBST(NM,I1,I2,J1,J2,K1,K2,CR_INDEX,OBST_INDEX)

USE MESH_POINTERS
USE GEOMETRY_FUNCTIONS, ONLY : BLOCK_CELL
INTEGER :: I1,I2,J1,J2,K1,K2,I,J,K
INTEGER, INTENT(IN) :: NM,CR_INDEX,OBST_INDEX
LOGICAL :: CREATE,REMOVE

CALL POINT_TO_MESH(NM)

! Indicate whether to create or remove the obstruction.

OBST_CREATED_OR_REMOVED = .TRUE.
REMOVE = .FALSE.
CREATE = .FALSE.
IF (CR_INDEX==0) REMOVE = .TRUE.
IF (CR_INDEX==1) CREATE = .TRUE.
IF (REMOVE) OBSTRUCTION(OBST_INDEX)%SCHEDULED_FOR_REMOVAL  = .TRUE.
IF (CREATE) OBSTRUCTION(OBST_INDEX)%SCHEDULED_FOR_CREATION = .TRUE.

! Blank or unblank cells that make up the OBSTruction

IF (I1/=I2 .AND. J1/=J2 .AND. K1/=K2) CALL BLOCK_CELL(NM,I1+1,I2,J1+1,J2,K1+1,K2,CR_INDEX,OBST_INDEX)

! If the OBSTruction is to be removed, set density and mass fractions to ambient value

IF (REMOVE) THEN
   DO K=K1+1,K2
      DO J=J1+1,J2
         DO I=I1+1,I2
            RHOS(I,J,K) = RHO_0(K)
            RHO(I,J,K)  = RHO_0(K)
            IF (SOLID_PHASE_ONLY) TMP(I,J,K) = TMP_0(K)
            ZZ(I,J,K,1:N_TRACKED_SPECIES)  = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
            ZZS(I,J,K,1:N_TRACKED_SPECIES) = SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
         ENDDO
      ENDDO
   ENDDO
ENDIF

END SUBROUTINE CREATE_OR_REMOVE_OBST


!> \brief Re-assign wall boundaries after create or removal of obstructions
!> \param T Current time (s)
!> \param NM Mesh number

SUBROUTINE REASSIGN_WALL_CELLS(T,NM)

USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY : CURRENT_TIME
INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: T
INTEGER :: I1,I2,J1,J2,K1,K2,I,J,K,IW,ICG,IC,OBST_INDEX,NOM,IIO,JJO,KKO
REAL(EB) :: TNOW
LOGICAL :: CREATE,REMOVE
TYPE (OBSTRUCTION_TYPE), POINTER :: OB
TYPE (WALL_TYPE), POINTER :: WC
TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (MESH_TYPE), POINTER :: MM

TNOW = CURRENT_TIME()

CALL POINT_TO_MESH(NM)

DO IW=1,N_EXTERNAL_WALL_CELLS
   EWC => MESHES(NM)%EXTERNAL_WALL(IW)
   NOM = EWC%NOM
   IF (NOM==0) CYCLE
   WC => MESHES(NM)%WALL(IW)
   BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)
   MM  => MESHES(NOM)
   IIO = EWC%IIO_MIN
   JJO = EWC%JJO_MIN
   KKO = EWC%KKO_MIN
   IF (WC%OBST_INDEX==0 .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY .AND. .NOT.MM%CELL(MM%CELL_INDEX(IIO,JJO,KKO))%SOLID) THEN
      IC  = CELL_INDEX(BC%II ,BC%JJ ,BC%KK )
      ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG)
      IF (CELL(ICG)%SOLID) CYCLE
      OBST_INDEX = 0
      REMOVE = .FALSE.
      CREATE = .FALSE.
      CALL GET_BOUNDARY_TYPE
   ENDIF
ENDDO

! Loop over all obstructions in the current mesh and initialize newly exposed or covered wall cell faces

OBSTRUCTION_LOOP: DO OBST_INDEX=1,N_OBST

OB => OBSTRUCTION(OBST_INDEX)

IF (.NOT.OB%SCHEDULED_FOR_REMOVAL .AND. .NOT.OB%SCHEDULED_FOR_CREATION) CYCLE OBSTRUCTION_LOOP

REMOVE = .FALSE. ; CREATE = .FALSE.
IF (OB%SCHEDULED_FOR_REMOVAL) THEN
   REMOVE = .TRUE.
   OB%SCHEDULED_FOR_REMOVAL = .FALSE.
ENDIF
IF (OB%SCHEDULED_FOR_CREATION) THEN
   CREATE = .TRUE.
   OB%SCHEDULED_FOR_CREATION = .FALSE.
ENDIF

I1 = OB%I1 ; I2 = OB%I2 ; J1 = OB%J1 ; J2 = OB%J2 ; K1 = OB%K1 ; K2 = OB%K2

DO K=K1+1,K2
   DO J=J1+1,J2
      IC  = CELL_INDEX(I1+1,J,K)
      ICG = CELL_INDEX(I1  ,J,K)
      IW  = CELL(ICG)%WALL_INDEX( 1)
      IF (IW>0 .AND. I1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I1  ,J,K)
      ICG = CELL_INDEX(I1+1,J,K)
      IW  = CELL(ICG)%WALL_INDEX(-1)
      IF (IW>0 .AND. I1<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2  ,J,K)
      ICG = CELL_INDEX(I2+1,J,K)
      IW  = CELL(ICG)%WALL_INDEX(-1)
      IF (IW>0 .AND. I2<IBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I2+1,J,K)
      ICG = CELL_INDEX(I2  ,J,K)
      IW  = CELL(ICG)%WALL_INDEX( 1)
      IF (IW>0 .AND. I2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the y boundaries of the OBSTruction

DO K=K1+1,K2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J1+1,K)
      ICG = CELL_INDEX(I,J1  ,K)
      IW  = CELL(ICG)%WALL_INDEX( 2)
      IF (IW>0 .AND. J1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J1  ,K)
      ICG = CELL_INDEX(I,J1+1,K)
      IW  = CELL(ICG)%WALL_INDEX(-2)
      IF (IW>0 .AND. J1<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2  ,K)
      ICG = CELL_INDEX(I,J2+1,K)
      IW  = CELL(ICG)%WALL_INDEX(-2)
      IF (IW>0 .AND. J2<JBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J2+1,K)
      ICG = CELL_INDEX(I,J2  ,K)
      IW  = CELL(ICG)%WALL_INDEX( 2)
      IF (IW>0 .AND. J2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Process the z boundaries of the OBSTruction

DO J=J1+1,J2
   DO I=I1+1,I2
      IC  = CELL_INDEX(I,J,K1+1)
      ICG = CELL_INDEX(I,J,K1  )
      IW  = CELL(ICG)%WALL_INDEX( 3)
      IF (IW>0 .AND. K1>0)    CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K1  )
      ICG = CELL_INDEX(I,J,K1+1)
      IW  = CELL(ICG)%WALL_INDEX(-3)
      IF (IW>0 .AND. K1<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2  )
      ICG = CELL_INDEX(I,J,K2+1)
      IW  = CELL(ICG)%WALL_INDEX(-3)
      IF (IW>0 .AND. K2<KBAR) CALL GET_BOUNDARY_TYPE
      IC  = CELL_INDEX(I,J,K2+1)
      ICG = CELL_INDEX(I,J,K2  )
      IW  = CELL(ICG)%WALL_INDEX( 3)
      IF (IW>0 .AND. K2>0)    CALL GET_BOUNDARY_TYPE
   ENDDO
ENDDO

! Nullify block edges on blockage that is to be removed

DO K=K1,K2
   DO J=J1,J2
      IF (J>J1) CALL REDEFINE_EDGE(I1,J,K,2)
      IF (J>J1) CALL REDEFINE_EDGE(I2,J,K,2)
      IF (K>K1) CALL REDEFINE_EDGE(I1,J,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I2,J,K,3)
   ENDDO
ENDDO

DO K=K1,K2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J1,K,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J2,K,1)
      IF (K>K1) CALL REDEFINE_EDGE(I,J1,K,3)
      IF (K>K1) CALL REDEFINE_EDGE(I,J2,K,3)
   ENDDO
ENDDO

DO J=J1,J2
   DO I=I1,I2
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K1,1)
      IF (I>I1) CALL REDEFINE_EDGE(I,J,K2,1)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K1,2)
      IF (J>J1) CALL REDEFINE_EDGE(I,J,K2,2)
   ENDDO
ENDDO

ENDDO OBSTRUCTION_LOOP

T_USED(6) = T_USED(6) + CURRENT_TIME() - TNOW
CONTAINS

!> \brief Determine the type and other properties of a newly exposed wall cell

SUBROUTINE GET_BOUNDARY_TYPE

INTEGER :: IOR,IIG,JJG,KKG,IW_OLD,IERR,PRESSURE_BC_TYPE,ICG_OLD,II
TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1,B1_OLD
TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D_OLD
TYPE (WALL_TYPE), POINTER :: WC_OLD
TYPE (EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE (STORAGE_TYPE), POINTER :: OS
TYPE (SURFACE_TYPE), POINTER :: SF

WC => MESHES(NM)%WALL(IW)
BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)

IF (IW<=N_EXTERNAL_WALL_CELLS .AND. REMOVE) THEN
   EWC => MESHES(NM)%EXTERNAL_WALL(IW)
   WC%BOUNDARY_TYPE = SOLID_BOUNDARY
   CELL(IC)%SOLID       = .TRUE.
   IF (EWC%SURF_INDEX_ORIG==MIRROR_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = MIRROR_BOUNDARY
      WC%SURF_INDEX    = MIRROR_SURF_INDEX
      CELL(IC)%SOLID   = .TRUE.
      RETURN
   ENDIF
   IF (EWC%SURF_INDEX_ORIG==OPEN_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = OPEN_BOUNDARY
      WC%SURF_INDEX    = OPEN_SURF_INDEX
      CELL(IC)%SOLID   = .FALSE.
   ENDIF
   IF (EWC%SURF_INDEX_ORIG==INTERPOLATED_SURF_INDEX) THEN
      WC%BOUNDARY_TYPE = INTERPOLATED_BOUNDARY
      WC%SURF_INDEX    = INTERPOLATED_SURF_INDEX
      CELL(IC)%SOLID   = .FALSE.
      RETURN
   ENDIF
ENDIF

IF (IW>N_EXTERNAL_WALL_CELLS) THEN
   IF (WC%OBST_INDEX>0 .AND.      OBSTRUCTION(WC%OBST_INDEX)%HIDDEN .AND. .NOT.CELL(IC)%SOLID ) WC%BOUNDARY_TYPE = NULL_BOUNDARY
   IF (WC%OBST_INDEX>0 .AND. .NOT.OBSTRUCTION(WC%OBST_INDEX)%HIDDEN .AND. .NOT.CELL(ICG)%SOLID) WC%BOUNDARY_TYPE = SOLID_BOUNDARY
   IF (CELL(ICG)%SOLID) WC%BOUNDARY_TYPE = NULL_BOUNDARY
ENDIF

IF (CREATE) THEN
   IF (CELL(ICG)%SOLID) THEN
      WC%BOUNDARY_TYPE = NULL_BOUNDARY
   ELSE
      WC%BOUNDARY_TYPE = SOLID_BOUNDARY
      B1 => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
      IF (B1%T_IGN<T) B1%T_IGN = T
   ENDIF
ENDIF

! Initialize solid properties included in the WALL derived type variable

IF (WC%BOUNDARY_TYPE/=NULL_BOUNDARY) THEN
   IF (IW<=N_EXTERNAL_WALL_CELLS) THEN
      WC%OBST_INDEX = 0
      EWC => EXTERNAL_WALL(IW)
      WC%SURF_INDEX = EWC%SURF_INDEX_ORIG
      PRESSURE_BC_TYPE = EWC%PRESSURE_BC_TYPE  ! Save this parameter and restore it after the call to INIT_WALL_CELL
   ENDIF
   IF (CELL(IC)%OBST_INDEX>0) THEN
      WC%OBST_INDEX = CELL(IC)%OBST_INDEX
      WC%SURF_INDEX = OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX(BC%IOR)
   ELSEIF (CREATE .AND. OBST_INDEX>0) THEN
      WC%OBST_INDEX = OBST_INDEX
      WC%SURF_INDEX = OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX(BC%IOR)
   ENDIF
   IF (OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX_INTERIOR>0) WC%SURF_INDEX = OBSTRUCTION(WC%OBST_INDEX)%SURF_INDEX_INTERIOR
   CALL INIT_WALL_CELL(NM,BC%II,BC%JJ,BC%KK,WC%OBST_INDEX,IW,BC%IOR,WC%SURF_INDEX,IERR,T)
   WC => MESHES(NM)%WALL(IW)
   IF (IW<=N_EXTERNAL_WALL_CELLS) EWC%PRESSURE_BC_TYPE = PRESSURE_BC_TYPE
ENDIF

! Special cases 1: BURNed_AWAY obstruction exposes a surface that also burns, in which case the surface is to ignite immediately.
! Special cases 2: HT3D solid shifts the position of the burned away surface to the exposed surface position.

SF => SURFACE(WC%SURF_INDEX)
IF (REMOVE .AND. ( (SF%THERMAL_BC_INDEX==THERMALLY_THICK.AND.(SF%VARIABLE_THICKNESS.OR.SF%HT_DIM>1)) &
                  .OR. SF%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED ) ) THEN
   BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)
   IIG = BC%IIG
   JJG = BC%JJG
   KKG = BC%KKG
   IOR = BC%IOR
   ICG_OLD  = 0
   SELECT CASE(IOR)
      CASE(-1) ; IF (IIG>1)    ICG_OLD = CELL_INDEX(IIG-1,JJG,KKG)
      CASE( 1) ; IF (IIG<IBAR) ICG_OLD = CELL_INDEX(IIG+1,JJG,KKG)
      CASE(-2) ; IF (JJG>1)    ICG_OLD = CELL_INDEX(IIG,JJG-1,KKG)
      CASE( 2) ; IF (JJG<JBAR) ICG_OLD = CELL_INDEX(IIG,JJG+1,KKG)
      CASE(-3) ; IF (KKG>1)    ICG_OLD = CELL_INDEX(IIG,JJG,KKG-1)
      CASE( 3) ; IF (KKG<KBAR) ICG_OLD = CELL_INDEX(IIG,JJG,KKG+1)
   END SELECT
   IW_OLD = CELL(ICG_OLD)%WALL_INDEX(-IOR)
   IF (IW_OLD>0) THEN
      WC_OLD => MESHES(NM)%WALL(IW_OLD)
      IF (SF%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN
         B1     => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
         B1_OLD => MESHES(NM)%BOUNDARY_PROP1(WC_OLD%B1_INDEX)
         IF (WC_OLD%SURF_INDEX==WC%SURF_INDEX) B1%T_IGN = B1_OLD%T_IGN
      ELSEIF (.NOT.CELL(ICG_OLD)%SOLID .AND. .NOT.CELL(ICG)%SOLID .AND. CELL(IC)%SOLID .AND. &
              SUM(BOUNDARY_ONE_D(WC_OLD%OD_INDEX)%N_LAYER_CELLS(:))>0) THEN
         WC%OD_INDEX = WC_OLD%OD_INDEX
         WC%BOUNDARY_TYPE = SOLID_BOUNDARY
         ONE_D_OLD => MESHES(NM)%BOUNDARY_ONE_D(WC_OLD%OD_INDEX)
         IF (ONE_D_OLD%BACK_MESH>0 .AND. ONE_D_OLD%BACK_MESH/=NM) THEN
            OS => OMESH(ONE_D_OLD%BACK_MESH)%WALL_SEND_BUFFER
            DO II=1,OS%N_ITEMS
               IF (OS%ITEM_INDEX(II)==IW_OLD) OS%ITEM_INDEX(II) = IW
            ENDDO
         ENDIF
      ENDIF
   ENDIF
ENDIF

END SUBROUTINE GET_BOUNDARY_TYPE


!> \brief Change a few properties of the EDGEs that have been exposed or covered up by a blockage
!> \param II x-index of edge
!> \param JJ y-index of edge
!> \param KK z-index of edge
!> \param IEC Edge index: 1=x, 2=y, 3=z

SUBROUTINE REDEFINE_EDGE(II,JJ,KK,IEC)

INTEGER :: IE,II,JJ,KK,IEC

SELECT CASE(IEC)
   CASE(1)
      IE = CELL(CELL_INDEX(II,JJ,KK))%EDGE_INDEX( 4)
   CASE(2)
      IE = CELL(CELL_INDEX(II,JJ,KK))%EDGE_INDEX( 8)
   CASE(3)
      IE = CELL(CELL_INDEX(II,JJ,KK))%EDGE_INDEX(12)
END SELECT

END SUBROUTINE REDEFINE_EDGE

END SUBROUTINE REASSIGN_WALL_CELLS


!> \brief Generate random noise at the start of the simulation
!> \param NM Mesh number

SUBROUTINE INITIAL_NOISE(NM)

USE MESH_POINTERS
REAL     :: RN2
REAL(EB) :: RN
INTEGER  :: I,J,K,SIZE_RND,IZERO
INTEGER, DIMENSION(:), ALLOCATABLE :: SEED_RND
INTEGER, INTENT(IN) :: NM

! Waste a few calls to RANDOM_NUMBER to avoid generating the exact same sequence on each mesh

CALL RANDOM_SEED(SIZE=SIZE_RND)
ALLOCATE(SEED_RND(SIZE_RND),STAT=IZERO)
CALL CHKMEMERR('INITIAL_NOISE','SEED_RND',IZERO)
SEED_RND = 2819 * 13*NM + RND_SEED
CALL RANDOM_SEED(PUT=SEED_RND)
DEALLOCATE(SEED_RND)

DO I=1,NM
   CALL RANDOM_NUMBER(RN2)
ENDDO

IF (.NOT. NOISE) RETURN

! Point to local mesh variables

CALL POINT_TO_MESH(NM)

! Add random vorticity to cells that are not bounding solid surfaces

DO K=1,KBM1
   DO J=1,JBM1
      DO I=1,IBAR
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDY(J)
         W(I,J+1,K) = W(I,J+1,K) + RN*RDY(J+1)
         V(I,J,K)   = V(I,J,K)   + RN*RDZ(K)
         V(I,J,K+1) = V(I,J,K+1) - RN*RDZ(K+1)
      ENDDO
   ENDDO
ENDDO
DO K=1,KBM1
   DO J=1,JBAR
      DO I=1,IBM1
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         W(I,J,K)   = W(I,J,K)   - RN*RDX(I)*R(I)*RRN(I)
         W(I+1,J,K) = W(I+1,J,K) + RN*RDX(I+1)*R(I)*RRN(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDZ(K)
         U(I,J,K+1) = U(I,J,K+1) - RN*RDZ(K+1)
      ENDDO
   ENDDO
ENDDO
DO K=1,KBAR
   DO J=1,JBM1
      DO I=1,IBM1
         CALL RANDOM_NUMBER(RN2)
         RN=REAL(RN2,EB)
         RN = NOISE_VELOCITY*(-1._EB + 2._EB*RN)*CELL_SIZE
         V(I,J,K)   = V(I,J,K)   - RN*RDX(I)
         V(I+1,J,K) = V(I+1,J,K) + RN*RDX(I+1)
         U(I,J,K)   = U(I,J,K)   + RN*RDY(J)
         U(I,J+1,K) = U(I,J+1,K) - RN*RDY(J+1)
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE INITIAL_NOISE


!> \brief Read UVW file
!> \param NM Mesh number
!> \param FN_UVW File name

SUBROUTINE UVW_INIT(NM,FN_UVW)

USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER,SHUTDOWN
INTEGER  :: I,J,K,II,JJ,KK,IW,IOR,LU_UVW,IERROR,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
INTEGER, INTENT(IN) :: NM
CHARACTER(80), INTENT(IN) :: FN_UVW
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
TYPE(WALL_TYPE), POINTER :: WC
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1

CALL POINT_TO_MESH(NM)

LU_UVW = GET_FILE_NUMBER()
OPEN(UNIT=LU_UVW,FILE=FN_UVW,FORM='FORMATTED',STATUS='OLD',IOSTAT=IERROR)

IF (IERROR/=0) THEN
   WRITE(MESSAGE,'(A,I0,A,A)') 'ERROR(439): MESH ',NM,', UVWFILE ',TRIM(FN_UVW),' does not exist.'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
   RETURN
ENDIF

IF (PERIODIC_TEST==2) THEN
   IMIN = 1
   IMAX = IBAR
   JMIN = 1
   JMAX = JBAR
   KMIN = 1
   KMAX = KBAR
ELSE
   READ(LU_UVW,*) IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
   IMIN = MAX(0,IMIN)
   IMAX = MIN(IBAR,IMAX)
   JMIN = MAX(0,JMIN)
   JMAX = MIN(JBAR,JMAX)
   KMIN = MAX(0,KMIN)
   KMAX = MIN(KBAR,KMAX)
ENDIF
DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         READ(LU_UVW,*,IOSTAT=IERROR) U(I,J,K),V(I,J,K),W(I,J,K)
         IF (IERROR/=0) THEN
            U(I,J,K)=0._EB
            V(I,J,K)=0._EB
            W(I,J,K)=0._EB
         ENDIF
      ENDDO
   ENDDO
ENDDO

CLOSE(LU_UVW)

IF (PERIODIC_TEST==2) THEN
   U(0,:,:) = U(IBAR,:,:)
   V(:,0,:) = V(:,JBAR,:)
   W(:,:,0) = W(:,:,KBAR)
ENDIF

US=U
VS=V
WS=W

! Set normal velocity on external and internal boundaries (follows divg)

DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
   WC => WALL(IW)
   BC => BOUNDARY_COORD(WC%BC_INDEX)
   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
   IOR = BC%IOR
   II  = BC%II
   JJ  = BC%JJ
   KK  = BC%KK
   SELECT CASE(IOR)
      CASE( 1) ; B1%U_NORMAL_S = -U(II,JJ,KK)
      CASE(-1) ; B1%U_NORMAL_S =  U(II-1,JJ,KK)
      CASE( 2) ; B1%U_NORMAL_S = -V(II,JJ,KK)
      CASE(-2) ; B1%U_NORMAL_S =  V(II,JJ-1,KK)
      CASE( 3) ; B1%U_NORMAL_S = -W(II,JJ,KK)
      CASE(-3) ; B1%U_NORMAL_S =  W(II,JJ,KK-1)
   END SELECT
   B1%U_NORMAL = B1%U_NORMAL_S
ENDDO

END SUBROUTINE UVW_INIT


!> \brief Read TMP file
!> \param NM Mesh number
!> \param FN_TMP File name

SUBROUTINE TMP_INIT(NM,FN_TMP)

USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER,SHUTDOWN
USE RADCONS, ONLY: UIIDIM
INTEGER  :: I,J,K,LU_TMP,IERROR,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
INTEGER, INTENT(IN) :: NM
CHARACTER(80), INTENT(IN) :: FN_TMP
CHARACTER(MESSAGE_LENGTH) :: MESSAGE

CALL POINT_TO_MESH(NM)

LU_TMP = GET_FILE_NUMBER()
OPEN(UNIT=LU_TMP,FILE=FN_TMP,FORM='FORMATTED',STATUS='OLD',IOSTAT=IERROR)

IF (IERROR/=0) THEN
   WRITE(MESSAGE,'(A,I0,3A)') 'ERROR(440): MESH ',NM,', TMPFILE ',TRIM(FN_TMP),' does not exist.'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
   RETURN
ENDIF

READ(LU_TMP,*) IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
IMIN = MAX(1,IMIN)
IMAX = MIN(IBAR,IMAX)
JMIN = MAX(1,JMIN)
JMAX = MIN(JBAR,JMAX)
KMIN = MAX(1,KMIN)
KMAX = MIN(KBAR,KMAX)

DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         READ(LU_TMP,*,IOSTAT=IERROR) TMP(I,J,K)
         IF (IERROR/=0) TMP(I,J,K)=0._EB
      ENDDO
   ENDDO
ENDDO

CLOSE(LU_TMP)

! update density field

DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         RHO(I,J,K)  = P_0(K)/(TMP(I,J,K)*RSUM(I,J,K))
         RHOS(I,J,K) = RHO(I,J,K)
         IF (RADIATION) THEN
            UII(I,J,K) = 4._EB*SIGMA*TMP(I,J,K)**4
            UIID(I,J,K,1:UIIDIM) = UII(I,J,K)/REAL(UIIDIM,EB)
         ENDIF
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE TMP_INIT 


!> \brief Read SPEC file
!> \param NM Mesh number
!> \param FN_SPEC File name

SUBROUTINE SPEC_INIT(NM,FN_SPEC)

USE MESH_POINTERS
USE COMP_FUNCTIONS, ONLY: GET_FILE_NUMBER,SHUTDOWN
USE RADCONS, ONLY: UIIDIM
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT,GET_REALIZABLE_MF
INTEGER  :: I,J,K,N,LU_SPEC,IERROR,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES)
INTEGER, INTENT(IN) :: NM
CHARACTER(80), INTENT(IN) :: FN_SPEC
CHARACTER(MESSAGE_LENGTH) :: MESSAGE

CALL POINT_TO_MESH(NM)

LU_SPEC = GET_FILE_NUMBER()
OPEN(UNIT=LU_SPEC,FILE=FN_SPEC,FORM='FORMATTED',STATUS='OLD',IOSTAT=IERROR)

IF (IERROR/=0) THEN
   WRITE(MESSAGE,'(A,I0,3A)') 'ERROR(441): MESH ',NM,', SPECFILE ',TRIM(FN_SPEC),' does not exist.'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
   RETURN
ENDIF

READ(LU_SPEC,*) IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
IMIN = MAX(1,IMIN)
IMAX = MIN(IBAR,IMAX)
JMIN = MAX(1,JMIN)
JMAX = MIN(JBAR,JMAX)
KMIN = MAX(1,KMIN)
KMAX = MIN(KBAR,KMAX)

DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         READ(LU_SPEC,*,IOSTAT=IERROR) ( ZZ(I,J,K,N), N=1,N_TRACKED_SPECIES )
         IF (IERROR/=0) ZZ(I,J,K,1:N_TRACKED_SPECIES)=0._EB
      ENDDO
   ENDDO
ENDDO

CLOSE(LU_SPEC)

! update density field

DO K=KMIN,KMAX
   DO J=JMIN,JMAX
      DO I=IMIN,IMAX
         ! Check realizability of input mass fractions
         ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)
         CALL GET_REALIZABLE_MF(ZZ_GET)
         ! Compute molecular weight term RSUM=R0*SUM(Y_i/M_i)
         ZZ(I,J,K,1:N_TRACKED_SPECIES) = ZZ_GET(1:N_TRACKED_SPECIES)
         CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(I,J,K))
         RHO(I,J,K)  = P_0(K)/(TMP(I,J,K)*RSUM(I,J,K))
         RHOS(I,J,K) = RHO(I,J,K)
         IF (RADIATION) THEN
            UII(I,J,K) = 4._EB*SIGMA*TMP(I,J,K)**4
            UIID(I,J,K,1:UIIDIM) = UII(I,J,K)/REAL(UIIDIM,EB)
         ENDIF
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE SPEC_INIT


END MODULE INIT
